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I .  INTRODUCTION 


The  objective  of  this  research  effort  is  to  develop  a  predictive 
numerical  2-D  model  which  meets  the  following  twenty-one  criteria: 


1.  Realistically  models  the  geometry  and  velocity  distribution  of 
any  stretching  rod  associated  with  state-of-the-art  self-forging 
fragments  (SSF)  and  shaped-charge  warheads. 

2.  Geometrically  models  any  zero-obliquity  target  containing  horizon¬ 
tal  or  vertical  layers,  multiple  materials,  free  surfaces  and 
possible  air  spaces. 

3.  Geometrically  formulated  so  that  the  stretching  rod  and  target 
model  can  be  readily  extended  to  oblique  impacts. 

4.  Includes  all  thermal  properties  of  any  rod  or  non-chemically 
reactive  target-material. 

5.  Equations-of-stace  which  account  for  phase  changes. 

6.  Compressible  and  viscous  effects  considered  for  fluids. 

7.  Elastic-viscoplastic  constitutive  equation  used  for  solid 
materials . 

8.  Strain-displacement  relationship  includes  large  rotations  and 
extensions . 

9.  Primary  shocks  are  included  in  the  model. 

10.  Time-dependent  fracture  criteria  is  used  for  solid  and  liquid 
materials . 

11.  Temperatures  within  the  solid  target  and  rod  can  be  determined 
if  initial  termperatures  and  convective  film  coefficients  are 
known. 

12.  Computation  times  should  be  less  than  10  minutes  for  a  2-D 
penetration  process  which  involves  500  ys  of  penetration  time. 

13.  Quantitative  agreement  with  existing  penetration  and  penetration 
vs.  time  data  for  shaped-charge  jets. 

14.  Quantitative  agreement  with  data  for  radial  target  response  for 
monolithic  or  layered  metallic  targets  impacted  by  a  stretching 
rod . 

15.  Quantitative  agreement  with  experimental  phenomenon  associated 
with  radial  target  influence  on  axial  penetration. 
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16.  Qualitatively  predicts  experimental  phenomenon  associated  with 
"particle"  jets. 

17.  Qualitatively  predicts  experimental  phenomenon  associated  with 
confined  columns. 
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18.  Qualitatively  predicts  experimental  phenomenon  associated  with 
ejected  material. 

19.  Qualitatively  predicts  experimental  phenomenon  associated  with 
heating  of  target  material  by  a  penetrating  jet. 

20.  Qualitatively  predicts  experimental  phenomenon  associated  with 
the  penetration  of  brittle  hard  materials. 

21.  Qualitatively  predicts  experimental  phenomena  associated  with 
the  influence  of  lateral  dimensions  and  confinement  upon  the 
behavior  of  brittle  hard  materials. 

This  2-D  model  is  to  be  developed  over  a  three-year  time-span,  and 
this  report  presents  the  results  of  the  effort  for  the  first  year. 

Obviously,  in  order  to  meet  the  above  criteria,  one  must  insure  that 
the  pertinent  physical  parameters  are  included  in  the  model,  and  that  the 
laws  of  physics  are  included  in  a  consistent  manner.  Therefore,  this 
work  concentrated  upon  the  establishment  of  a  self-consistent  set  of 
equations.  However,  independent  and  University  sponsored  research  was 
concerned  with  the  review  of  state-of-the-art  solution  techniques  used 
in  existing  structural  and  hydrodynamic  codes.  Since  this  independent 
research  had  a  bearing  on  the  development  of  the  equations,  the  review 
will  be  discussed  first. 


II.  OVERALL  SOLUTION  SCHEME 

Recent  publications  involving  non-linear  dynamics  of  solids,^  the 

2.  3  * 

flow  of  solids  during  extrusion,  numerical  heat-transfer,  and  computational 


^■J.  T.  Oden  and  G. 
:-^ECHA::ICS,  volume 
.Vev  Jersey,  1534. 


F.  Carey,  FINITE  ELEMENTS:  SPECIAL  PROBLEMS  IN  SOLID 
V  of  the  Texas  Finite  Element  Series,  Prenviae-Hall  Inc. , 


"G.  C.  Zienkiewiaz  and  P.  N.  Godbole,  "Flow  of  Plastic  and  Visco-Plastic 
Solids  with  Special  Reference  to  Extrusion  and  Forming  Processes ,  "  Int.  J. 
Nian.  Mezh.  Engr. ,  S,  3-16,  1574. 

^ - 

^0.  C.  Zienkiewicz  and  C.  J.  Parekh,  "Transient  Field  Problem.s:  pjo- 


jzmeyiszcna.-  ay 
Elem.enzs ,  "  In: 


Three- Gim.ensional  Analysis  by  Isoparamezrzc  Fzmte 


IJumerz^ 


Methods  Engr. , 


1570. 


fluid-mechanics,**  show  chat  the  finite-element  technique  (mathematically 
known  as  the  Method  of  Weighted  Residuals)  can  yield  good  approximations. 
Most  importantly,  a  "good"  approximate  solution  can  be  obtained  using  a 
relatively  small  number  of  judiciously  selected  "super  elements."  Also, 
the  problem  appears  to  be  well-posed  if  velocities  and/or  temperatures  are 
taken  as  the  "essential"  boundary  conditions,  with  the  stresses,  pressure 
and  heat  fluxes  taken  as  the  "natural"  boundary  conditions.  Hence,  these 
ideas  must  be  kept  in  mind  during  the  establishment  of  the  self-consistent 
set  of  equations. 

Figure  1  shows  an  assumed  system  associated  with  one  particular  region 
of  material.  This  region  of  material  is  similar  to  a  moving  control-volume 
with  a  time-dependent  shape  and  volume.  The  volume  and  shape  changes  due 
to  both  boundary  and  internal  conditions,  and  possible  mass  transfer  across 
its  boundaries.  Note  that  the  applied  boundary  conditions  can  involve 
velocity,  and  surface  pressures  (negative  tractions)  due  to  detonating 
explosives.  Figure  2  shows  a  schematic  diagram  of  a  region  with  a  sub- 
region  undergoing  a  phase-transition  due  to  the  applied  boundary  conditions. 
This  phase-transition  could  be  either  a  thermodynamic  transition,  e.g., 
solid  to  liquid  via  a  melting  region,^  or  a  polymorphic  transition  involving 
the  atomic  structure  of  the  solid  material.^ 

Since  the  thermal-mechanical-failure  properties  of  the  material  undergo 
drastic  changes  due  to  either  type  of  phase-transition,  the  author  chose  to 
represent  each  phase-region  of  material  as  a  "super-element"  of  the  finite- 
element  approximation.  Furthermore,  any  number  of  interacting  super¬ 
elements  can  be  used  to  represent  any  particular  problem  of  interest,  i.e., 
a  shaped -c;Tar.;e  jet  impacting  a  target,  or  a  shaped-charge  jet  being 
impacted  by  anotb.er  body. 

Note  that  this  usage  of  separate  elements  for  each  "phase”  region  can 
cause  the  "creation"  and  "annihilation"  of  elements  during  the  penetration 
orocess.  Mathematically ,  this  is  included  in  all  the  conservation  equations 
via  mass,  momentum  and  energy  transfer  across  the  phase-boundary  Spg  joining 
two  different  phase  regions.  Furthermore,  prior  to  a  phase-transition, 
certain  phase  regions  will  belong  to  a  "null-set,"  i.e.,  the  physical  size, 
mass,  momentum  ana  energy  are  identically  zero.  Numerically,  because  of  the 
finite  time-steps,  the  initial  zone  of  a  new  phase  is  mapped  from  the 
solution  at  the  previous  time  step,  e.g.,  the  geometric  locations  where 
T^Tm.  Because  of  certain  discontinuities  in  volume,  and  endothermic  or 
exothermic  entropy  c’nanges.  a  "refined"  phase-zone  may  be  required  via  an 
iterative  solution  at  that  particular  time-step. 


J .  Bake.',  'FINITE  ELEMENT  'lOMPi'IATIJNAL  FLUID  MECHANICS ,  Herr.is 


FublishiKa  Carp.,  New  Iovk, 


^L.  H.  VanVlack,  MA.IERIALS  FDR  ENGINEERING,  Chapters  5  and  6,  Add 
Wesley  Fublishina  Co.,  Heading,  Mass.,  1382. 
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The  above  process  of  creation  and  annihilation  of  elements  (hereafter 
referred  to  as  the  "C  and  A"  process)  is  somewhat  similar  to  the  process 
of  "rezoning"  used  for  highly  distorted  regions  in  Lagrangian  hydrocodes 
such  as  HEMP.  However,  in  the  C  and  A  process,  the  mass,  momentum  and 
energy  of  the  annihilated  material  is  mapped  directly  (within  the  accuracy 
of  the  finite-element  method)  into  the  created  material.  Conversely,  (to 
the  author's  knowledge)  the  mass,  momentum  and  energy  of  the  highly 
distorted  regions  are  "lost"  during  a  rezoning  process. 

It  appears  that  the  above  C  and  A  process  using  super-elements  differs 
considerably  from  current  finite-element  (hereafter  denoted  as  FE)  methods 
concerned  with  contact  and  indentation  problems  (see  review  in  reference  1) . 
Consequently,  a  few  comments  will  be  made  concerning  current  FE  methodology 
and  the  proposed  extension  to  the  C  and  A  process. 

The  most  commonly  used  technique  for  solids  involves  the  displacement 
formulation  which  requires  a  semi-inverse  method  where  the  functional  form 
of  the  displacement  solution,  within  the  FE  region  is  assumed  a  priori. 

However  in  penetration  problems,  initial  and/or  interface  velocities  are 
specified,  certain  velocities  are  knoi^n  (usually  zero  at  certain  boundaries), 
fluids  can  occur  within  the  target  and  in  many  problems  velocities  are  the 
desired  quantities.  Therefore,  the  author  chose  to  use  a  velocity  formula¬ 
tion*  for  the  FE  model.  Because  of  the  large  amount  of  existing  radio- 
graphic  and  framing-camera  data  and  various  hydrodynamic  code  calculations, 

I  believe  that  good  initial  guesses  on  velocity  distributions  can  be  made. 
Obviously,  if  a  situation  arises  where  no  information  exists  or  the  initial 
guesses  prove  to  be  "poor,"  any  material  region  can  be  subdivided  into 
numerous  interconnected  subdomains  i.e.,  more  than  one  super-element.  If 
the  solution  technique  is  properly  formulated,  one  should  be  able  to 
"converge"  to  a  solution  of  the  desired  accuracy  with  respect  to  a  pre¬ 
selected  measure  of  error.  Hence  the  solution  technique  must  contain  a 
consistent  "convergence"  criteria,  and  this  is  currently  being  studied  via 
independent  and  University-funded  research. 

As  will  be  pointed  out  later,  the  Method  of  Weighted  Residuals,  as 
implemented  by  the  FE  Method,  involves  the  solution  of  a  simultaneous  set 
of  non-linear  algebraic  equations,  whose  degree  of  non-linearity  depends 
upon  the  magnitude  of  the  deformations  and  the  constitutive  equation  for  the 
material.  The  review  indicated  that  the  most  popular  method  of  solving  this 
set  of  equations  is  the  following  approach: 

a.  drop  all  non-linear  terms  and  solve  the  simultaneous  linear  set 

b.  add  the  first  non-linear  term  to  the  set  of  equations,  in  the 
previous  linear  solution  and  then  iterate  until  a  certain 
criteria  of  convergence  is  achieved 

•k 

Displacements  are  kinematically  related  to  the  velocities  via  a  simple 
time-integration  over  the  velocity  history. 


c.  add  the  next  non-linear  term  to  the  set  of  equations  and  repeat 
step  b 

d.  continue  the  above  approach  until  all  non-linear  terms  have 
been  accounted  for. 

The  advantage  of  the  above  approach  is  that  the  addition  of  the  non¬ 
linear  terms  can  be  controlled  by  the  user.  Hence,  numerical  experiments 
can  be  conducted  to  ascertain  the  importance  of  the  various  non-linear 
contributions  and  those  deemed  of  lower-order  contributions  can  be  "turned- 
off"  by  the  user.  Therefore,  the  user  can  do  a  trade-off  between  accuracy 
of  solution  and  the  amount  of  computer  time  it  takes  to  achieve  an  approxi¬ 
mate  solution.  It  appears  that  nodal  velocities  are  an  appropriate  measure 
of  error  for  usage  in  the  above  iteration-scheme,  and  this  will  be  investi¬ 
gated  during  the  second  year  of  research. 

A  disadvantage  of  the  above  FE  method  is  that  the  iterative  solution 
of  simultaneous  equations  is  a  slow  process  when  there  is  a  large  number  of 
unknowns  in  the  so-called  bandwidth  or  wavefront.  Furthermore,  numerical 
fluid -mechanics  studies  (see  review  in  reference  4)  have  shown  that,  for  the 
same  number  of  unknowns,  it  takes  longer  to  solve  a  system  of  super-elements 
(so-called  quadratic  or  cubic  elements)  than  a  larger  system  of  linear 
elements  (such  as  used  in  the  EPIC  and  DYNA3D  programs).  Conversely, 
application  of  a  few  super-elements  to  solid  mechanics  problems^’®  have 
produced  "quicker"  solutions  when  the  problem  is  properly  formulated.**’® 
Therefore,  a  system  consisting  of  a  few  judiciously  selected  super-elements 
st\ould  provide  "quick"  approximate  solutions. 

Historically,  numerical  solutions  involving  shaped-charge  jet  penetra¬ 
tion  necessitate  a  large  number  of  zones  or  elements  because  of  the  small¬ 
ness  of  the  jet  (order  of  ram)  compared  with  the  target  thickness  (hundreds 
of  mm).  Furthermore,  very  large  velocity-gradients  occur  in  a  small  region 
(order  of  the  penetrator-diameter)  around  the  moving  penetrator/target 
interface.  This  poses  a  dilemma  in  that  our  "few  judiciously  selected  super¬ 
elements"  cannot  a  priori  model  this  condition. 

It  appears  to  the  author  that  the  following  approach  depicted  in 
Figures  3  through  7  offers  a  way  around  both  the  above  dilemma,  and  for 
implementing  the  C  and  A  process.  First,  we  recognize  that  the  vast 
majority  of  most  targets  undergo  linear  behavior  in  lateral  regions  which 
are  on  the  order  of  several  hole  diameters  away  from  the  path  of  the 
penetrator.  Hence,  these  outer  regions  could  be  modeled  using  a  few  super- 


'0.  C.  Zienkieuiaz,  THE  FINITE  ELEMENT  METHOD  IN  ENGINEERINQ  SCIENCE, 
Chapter  9,  MaGraw-Hill,  London,  1971. 

J.  Segerlind,  APPLIED  FINITE  ELEMENT  ANALYSIS,  Chapters  14-16,  Cohn 
Wiley  <3  Sons,  Inc.,  New  York,  1976. 

^C.  F.  Carey  and  C.  T.  Qden,  FINITE  ELEMENTS:  A  SECOND  COURSE,  Volume  II 
in  the  Texas  Finite  Element  Sertes,  Chapters  1  and  9,  Prenttae-Hall ,  Inc., 
New  Jersey,  1933. 
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FIGURE  4  Schematic  Overview  of  Deformations 
at  Three  Different  Times 


FIGURE  6  Active  Elements  at  Time 

Showing  the  Interface  Vector-Functions  g(t) 
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elements  on  either  side  of  the  penetrator-path.  Furthermore,  the  algebraic 
equations  are  essentially  linear  for  the  FE  unknowns  within  these  regions. 
Secondly,  the  "inner"  region  along  the  path  of  the  penetrator  is  the  only 
region  for  the  highly  non-linear  C  and  A  process.  Hence,  the  non-linear 
formulation  is  applied  only  to  this  "inner"  region.  Thirdly,  the  inner 
and  outer  regions  communicate  only  via  the  matching  of  velocities, 
tractions,  temperature  and  heat-fluxes  along  the  common  interface  (see 
Figure  3) .  This  matching  must  be  done  in  an  iterative  fashion,  but  does 
not  have  to  be  done  at  each  time-step.  Fourthly,  the  C  and  A  process 
progresses  along  the  penetrator-path  in  some  fashion  with  the  penetrator/ 
target  interface,  and  downstream  material  does  not  sense  the  penetration 
process  until  an  elastic  or  plastic  stress-wave  arrives  (see  Figure  4). 
Figures  5  through  7  illustrate  how  inner  super-elements  become  active,  i.e. , 
must  be  included  in  the  simultaneous  iterative  solution,  whereas  other 
inner  super-elements  become  inactive, with  the  motion  of  the  penetrator/ 
target  interface.  Note  that  this  active/ inactive  process  is  exactly 
analogous  to  Che  C  and  A  process,  and  all  iterative  interfaces  can  be 
analyzed  by  the  same  method. 

Numerically,  the  above  C  and  A  process  appears  to  offer  several 
advantages  over  current  FE  methodologies.  First,  a  priori  knowledge  of 
the  initial  location  and  material  associated  with  each  phantom  super- 
element  allows  the  user  to  "turn-on"  the  appropriate  non-linearities 
within  each  element,  which  in  general  can  differ  from  element  to  element. 
This  could  permit  a  substantial  savings  in  computational  ^time. 

Secondly,  it  is  numerically  faster  to  solve  N  sets  of  M  number  of 
simultaneous  equations  associated  with  each  super-element,  than  the  N*M 
simultaneous  equations  associated  with  all  N  super-elements.  Since  an 
iterative  solution  must  be  obtained  regardless  of  the  number  of  simul¬ 
taneous  equations,  minimal  solution  times  should  exist  for  certain  choices 
of  N  and  M.  Therefore,  a  trade-off  should  exist  between  the  number  of 
preselected  phantom  super-elements  and  the  number  of  nodes  associated  with 
each  super-element. 

Thirdly,  since  the  outer  super-elements  and  each  inner  super-element 
is  analyzed  or  "processed"  separately,  considerably  "faster"  solutions  are 
possible  by  using  computers  with  multiprocessors.  For  example,  a  co¬ 
processor  or  a  parallel  processor  could  be  used  to  run  the  inner  and  outer 
regions,  which  "calk  to  each  other"  after  a  preselecced  number  of 
iterations  or  time-steps.  This  may  be  extremely  important  in  the  near 
future  since  the  next  generation  of  "super  computers"  are  purported  to  be 
based  upon  a  large  number  of  interconnected  processors. 
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III.  DEVELOPMENT  OF  THE  PHYSICAL  MODEL 


Now  that  the  flavor  of  the  t-olution  scheme  has  been  presented,  we 
can  progress  to  the  assumptions  associated  with  the  selected  physical  model. 
The  corresponding  self-consistent  set  of  equations  must  include  the 
conservation  laws  of  physics  and  chemistry,  the  basic  principals  of 
thermodynamics  and  materials  science,  and  the  basic  postulates  of  con¬ 
tinuum  mechanics.  Presently,  the  model  will  assume  no  chemical  reactions 
and  hence  chemical  considerations  will  be  neglected.  Also,  the  model 
neglects  all  electromagnetic  forces  and  energy,  contains  no  failure  or 
fracture  criteria,  nor  any  consideration  of  solution  bifurcations,  i.e., 
buckling  phenomena  is  not  considered  explicitly. 

As  indicated  in  Figure  1,  both  the  global  inertial  coordinate  and  the 
attached  local,  fixed,  coordinate  systems  are  chosen  to  be  rectangular 
Cartesian  coordinates.  Hence,  all  the  equations  will  be  written  and  derived 
with  respect  to  these  coordinates.  Application  to  other  types  of  coordin¬ 
ate  systems  can  be  done  using  the  rules  of  calculus  transformations.  Note 
that,  while  the  local  coordinate  system  is  attached  to  some  material  point 
in  the  system  which  moves  globally  according  to  the  position  vector  R(t), 
the  direction  of  the  axes  always  corresponds  to  the  direction  of  the  global 
axes . 

Also  attached  to  the  material  system  is  a  curvilinear  material 

(so-called  Lagrangian)  coordinate  system  initially  parallel  to  the 
rectangular  Cartesian  system,  but  which  deforms  with  the  material.  Fur¬ 
thermore,  each  super-element  has  its  own  local  coordinate  system  used  for 
calculations.  After  each  calculation,  the  actual  location  of  the  total 
material  svstem  can  be  obtained  bv  mapping  back  to  the  inertial 
coordinate' system ,  and  this  yields  the  updated  nodal  locations.  It 

is  this  mapping  to  the  actual  location  chat  permits  one  to  observe  any 
possible  penetrator/ target  interactions  away  from  the  frontal  interface. 

Using  these  coordinates  and  including  the  possibility  of  phase 
transitions,  conservation  of  mass  yields  the  equation  presented  in 
Appendix  A.  This  equation  is  satisfied  (within  the  accuracy  of  the  mesh 
or  grid)  by  using  the  updated  Lagrangian  coordinates  and  the  previously 
discussed  C  and  A  approach. 

In  order  to  derive  the  remaining  equations,  we  must  decide  on  the 
magnitude  of  the  deformations  that  occur  between  time-steps.  The  usual 
infinitesimal  form  of  the  conservation  equations  requires  extremely 
small  time-steps  for  a  penetration  problem,  i.e.,  on  the  order  of 
-1  -2 

10  to  10  us.  This  would  require  a  large  amount  of  CPU  time  for  a  FE 
solution  to  a  realistic  problem.  Therefore,  since  relatively  large  time- 
steps  (order  of  usee)  are  desirable,  finite  deformations  are  quite 
possible  between  time-steps.  Hence,  large  deformation  theory''^  will  be 
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utilized  in  the  current  model  and  introduces  geometric  sources  of  non¬ 
linearity.  However,  the  advantage  of  large  deformation  theory  is  that 
the  geometric  deformations  do  not  force  one  to  use  small  time-steps  in  the 
solution.  Rather,  the  size  of  the  time  steps  is  dictated  by  either: 
a.  the  time  variation  of  the  boundary  conditions,  or  b.  the  user's 
interest  in  short  time  behavior,  i.e.,  "shock  propagation." 

Several  assumptions  must  be  introduced  into  the  large  deformation 
theory  in  order  to  simplify  the  elemental  form  of  the  conservation 
equations.  These  assumptions  involve  the  size  of  the  extensions  and 
shears  that  occur  between  the  time  steps.  This  size  is  assumed  to  be 
finite  but  small  compared  to  unity  and  one  radian  respectively.  Note 
chat  Che  theory  considers  material  rotations  which  could  be  larger  chan 
either  the  extensions  or  Che  shears.  Previous  hydrocode  and  experimental 
observations  indicate  that  a  material  "particle"  can  undergo  considerably 
large  rotations  at  the  target/penetrator  interface  region,  and  hence  no 
limitations*  were  applied  Co  the  material  rotations.  The  author  believes 
chat  Che  extension  and  shear  assumptions  are  not  severely  limiting 
because,  for  a  stretching  rod,  these  assumptions  apply  as  long  as  the 
following  two  conditions**  are  satisfied: 


Small  extensions 


Small  shears 


Ll  +  it  ] 


At  2 


where  At  =  time-step  size, 

c  =  time  of  the  calculation  as  measured  from  the  jet-formation  time. 

OB 

Using  Che  above  assumptions,  Che  linear  momentum  equation  takes  the 
form  given  in  Appendix  B.  The  body  force  terms  are  written  in  the  expanded 
form  since  thia  permits  a  visual  interpretation  of  the  various  contrib¬ 
utions,  and  is  easier  to  write  in  FORTRAN.  Also,  the  inertial  (body) 
forces  includes  all  forms  associated  with  rotational  and  translational 
motion  of  the  local  coordinate  system. 

Note  that,  if  the  body  forces  or  surface  tractions  are  either  uniform 
or  symmetrical  with  respect  to  the  local  coordinate  system,  a  geometric 
svmimetry  with  respect  to  the  same  axis  will  yield  a  zero  contribution  for 

*Ro tat  ions  a)  of  the  vectors  are  considered  during  time-steps,  but  their 
magnitudes  must  be  such  that  tanLw(c  +  At)  -  lij(t)]  =  [wCt  +  At)  -  ^(t)]. 

**  These  conditions  stem  from  the  assumption  that  the  radial  displacement 
rate  at  the  interface  is  of  the  same  order  as  the  rate  of  penetration. 

Also  it  is  assumed  chat  the  Cauchy  stress  is  approximately  equal  to  the 
first  Piola-Kirchhof f  stress. 
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the  corresponding  volume  or  surface  integral.  Also,  many  of  the 
volumetric  integrals  will  yield  geometric  terms,  =uch  as  the  polar 
moments  of  inertia,  when  the  other  integrands  are  constant  over  the 
volume . 

Appendix  C  gives  the  development  of  the  three  equations  associated 
with  the  conservation  of  angular  momentum.  Note  that,  whereas  the  point- 
wise  form  of  the  angular  momentum  equations  are  identically  satisfied  when¬ 
ever  the  linear  momentum  equations  are  satisfied,  the  weighted  integral  form 
requires  additional  symmetry  conditions.  Therefore,  these  angular  momentum 
equations  could  be  satisfied  by  either  assuming  that  these  symmetry 
conditions  are  fulfilled  by  the  deformation  conditions,  or,  by  solving  for 
the  additional  unknowns. 


The  last  conservation  equation  involves  energy  and  is  developed  in 
Appendix  D.  Because  of  the  assumed  general  behavior,  this  equation  is 
quite  complex  even  though  only  one  new  variable  (temperature)  is  introduced. 
Furthermore,  since  stored  thermal  energy  is  required,  some  form  of  thermo¬ 
dynamic  equation-of-state  (hereafter  denoted  as  EOS)  must  be  selected. 

Since  temperature  is  the  critical  variable  in  the  C  and  A  process,  all 
thermal-mechanical  properties  will  be  assumed  to  be  known  functions  of 
temperature  T.  Refereuccs  11-13  give  some  typical  types  of  analytical 
functions  used  for  these  properties. 


Several  new  terms  which  appear  in  the  energy  equation  are  the  strains 
the  heat  sources  Q  per  unit  mass,  and  tae  volume  fractions  fp  and  fg. 

Because  of  the  possibility  of  large  deformations,  tiie  itrain-displacement 
(Ui,  Uj ,  Uk.)  relationship  is 


,  3U.  3U.  3U,  3U 

c  =  i  (—i  +  __J.  +  X 

^ij  2  ^3x.  3x,  X.  3x  ^ 

J  ^  ^  j 


(1) 


where  i,  j,  k  =  1,  2,  3  and  repeated  indices  indicates  summation.  The  heat 
source  term  is  associated  with  the  latent  heats  corresponding  to  thermo¬ 
dynamic  and  polymorphic  phase-changes  within  the  material.  Most  of  these 
latent  heats  can  be  obtained  from  solid-state  reference  books.  The  volume 


*5.  K.  Godunov,  A.  F.  Demchuk,  N.  S.  Kozin  and  V.  I.  Mali,  "Intervclazion 
Formulas  for  Maxioell  Viscosity  of  Cervain  Metals  as  a  Function  of  Skear- 
Strain  Intensity  and  Temperature,  ”  translated  from  Zhumal  Frikl^dno  i 
Prikladnoi  Mekhaniki  i  Teknicheskoi  Fiziki  11974 j ,  ? ) enum  Pao lisnl n: 
Carp.,  1976. 
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■‘“J.  S.  Duvall  and  D.  P.  Dandekar,  "Theory  of  Equations  of  State:  Elastic- 
Plastic  Effects  II,"  USA  Ballistic  Research  Laboratories  Contract  Rorort 
Go.  106,  May  1975. 

^'fi.  3.  "alien,  T EERMODER AMISS ,  .Appendices  D  and  E,  -John  Wiley  and  Sons, 
Inc.  ,  lle'J  Tcrk,  4th  Printing ,  1963. 


fractions  fp  and  fg  refer  to  the  percentage  of  "plastic"  work  that  is 
irreversible  and  reversible,  with  many  processes  yielding  fe  of  only  a  few 
percent.*  However,  this  small  percentage  may  be  important  in  some  problems. 


The  two  remaining  items  that  must  be  selected  are  the  specific  forms 
of  the  EOS,  and  the  constitutive  equation  that  relates  stress  to  the 

other  variables. 


A.  Thermodynamic  Equations  of  State 

It  appears  that  most  hydrocodes  use  an  explicit  volume-dependent  form 
of  the  solid  EOS,  such  as  LASL's  equation^‘* 


P(p)  =  '^’s^  (Internal  Energy  -  1^^) 


(2) 


where 

P^(p)  =  pC“X/{ 1  -  X(S-l)}^, 

U_  =  Shock  Velocity  =  C  +  SU  , 
S  _  p 

Up  =  particle  velocity  =  v, 

X  =  (p/p,^)  - 


Y,  =  Gruneisen  coefficient. 


I.,  =  i  / p — ,  -  '  ■>  ,  C  =  Intercept  of  the  U  versus  U  data-line, 

n  1  \  (p/oo)-SX>  3  p 


p  =  density  in  the  deformed  state,  P(p)  =  pressure, 

p  =  reference  density,  S  =  slope  of  the  U  versus  U„  data-line. 

0  s  p 

Conversely,  many  of  tne  elastic-plastic  codes  use  a  Mie-Gruneisen  or 
Cristescu  form^^  which  explicitly  includes  the  influence  of  temperature, 
i  .e .  , 

PU,T)  =  P^(p)  +  Hq(p,T)  +  H^(p)T‘  (3) 


?iS‘ 


zbcrci 


tor'j  y  F!epcrt 


plasticit 
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used  to  represent 


where  Hq(p,T)  and  H^(p)  depend  upon  the  atomic  model^^"^^ 

Che  solid  state. 

Our  EOS  must  include  other  phases,  and  presently  there  appears  to  be 
two  choices  for  multi-thermodynamic  phase  EOS;  the  Tillocson^'  interpolation 
equation,  or  one  of  the  various  forms  of  the  GRAY^^  equation.  Tillotson's 
equation  is  a  two-phase  equation  based  upon  a  Thomas-Fermi-Dirac  atomic  model 
(five  parameters)  for  a  solid  (condensed)  state,  and  the  vaporized  (isentro- 
pically  expanded)  material  is  represented  by  a  two-parameter  equation  which 
interpolates  between  the  solid  state  and  an  ideal  gas.  GRAY  is  a  three- 
phase  equation  for  metals  based  upon  a  variety  of  physical  models  and  inter¬ 
polation  equations  (see  review  in  reference  19) . 

Neither  the  GRAY  nor  Che  Tillotson  equations  contain  polymorphic  phase- 
changes.*  (But  Che  author  believes  that  these  changes  can  be  included  by 
using  a  dual  Hugoniot  for  Pjj(p),  with  a  possible  volume  discontinuitv 
between  Che  dual  Hugoniot  relations.)  Cristescu  discusses  this  behavior  in 
reference  15,  and  some  specific  examples  will  be  considered  during  the 
research  effort  of  the  second  year. 

Since  BRL  personnel  have  used  the  GRAY  equation,  and  have  already 
generated  the  critical  thermodynamic  constants  associated  with  the  liquid- 
vapor,  mixed-phase  region  for  seventeen  materials, the  author  selected  the 
GRAY  form  of  Che  EOS.  After  reviewing  the  derivations  contained  in  reference 
19,  it  appears  chat  there  are  two  "implicit"  assumptions  contained  within 
the  equations  for  the  solid  region.  These  assumptions  are:  a.  the  initial 
temperature  of  the  solid  is  considerably  above  the  Debye  temperature  Tq,  and 
b.  the  energy  depends  only  on  the  thermal  and  volumetric  stored  energies. 
Since  Tq  could  be  well  above  room  temperature  (see  discussion  in  reference 
11)  Che  first  assumption  is  removed  from  the  solid  EOS  by  replacing  the  "3R  " 


■J .  ie  Bearjnent  2nd  J.  Leygonie ,  ’’Vapovizazion  of  'Jraniun  afzsr  Shock 
Loading,”  Proceedings  of  the  Fifth  Symposiwr’  (Inzematicnal)  on  Detonation, 
Pasadena,  California,  ACP-184-,  pp  547-553,  August  13-21,  1370. 
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J.  H.  Tillotson,  "Metal  Equations  of  State  for  Hypervelocity  Inpact," 
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the  GRAY  Equation  of  State,”  J.S.  ARPADCCt-I  Ballistic  Research  Laocratory 
Technical  .Report  Ax.3RL-TR-j2253 ,  August  133G. 
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The  HEMP  USER'S  MANU.YL,  Lawrence  Livermore  Laboratory  report  UCRL-51079 
Rev.  1,  Dec  1973,  indicates  a  "five-part  multiphase  equation  of  state," 


which  appears  to  be  for  pi.'ssibie  polymorphic  phase-changes  since  Che 
"phases”  are  determined  or.lv  by  the  range  of  the  densitv  change  X. 


terms  by  the  more  general  expression  f  c,, (DdT/W. 

o  ^ 

where , 

c^(T)  =  Temperature-dependent  coefficient  of  specific  heat  for 

R'  =  Universal  Gas  Constant  R/Molecular  weight 
Limit  c  (T)  3  3R. 

T  -  Td 

The  second  assumption  can  be  avoided  by  not  using  the  Gruneisen- like 
equation  with  parameter  r(T,V),*  but  rather  by  using  the  Hugoniot  defini¬ 
tions  (Equations  5,  6,  7  and  9  in  reference  19)  with  the  above  modified 
solid  EOS.  Therefore,  the  solid  EOS  takes  the  following  form, 

o  A 

P^(T,p)  =  ?^{p)  [1-  Yg(V)X/2]  +  Yg(V)p  Eq^(p) 


+  r.OTo  /  ogV 

3  2  e 

o  A 


where 


Eq^  =  assumed  by  the  present  author  to  be  the  linear  elastic 

T 

bulk-behavior  at  T^,  i.e.,  2K:(Tp)  /  a (T)dT/p  (T^) , 

< (Tq)  =  elastic  bulk-modulus  coefficient, 
o[(T)  =  unidirectional  coefficient  of  thermal  expansion, 

EqiCp)  =  Royce's  approximation-function 

I^X)  n  +  fx  +  (1-  ^)x2)  +  E.  (1+X„X), 


Aq  =  Lattice  constant, 

■  =0H  - 

(reference  19  uses  -Tp  {  3*8 . 134*10“^+150G^/W,5^} ) , 


It  appears  that  this  term  is  never  used  in  reference  19,  but  merely 
represents  a  formal  correction  term  to  a  Griineisen-like  EOS  for  a  solid. 


G’  =  atomic-weight  scaled  electronic  energy  coefficient, 

=  Ng  K  R'/2Ef, 
ge  =  Wa^' . 

K  =  Boltzmann's  constant, 

Ng  =  number  of  free  electrons, 

Ej  =  Fermi  Energy, 

=  26  (Ng 

Yg(V)  =  linear  volume-dependent  Griineisen  coefficient, 

Yg  =  electronic  gamma  (reference  19  assumes  y  =2/3 
whereas  references  15  and  16  assume  Yg  =  4) ■ 

With  the  above  modified  form  of  the  solid  EOS,  we  can  use  the  remaining 
GRAY  phase-equations, 

a.  Two-Phase  Melt  Region 

E^  (T,p)  =  E^(T,p)  +  v{T-  v(AT/2)}  (AS ' -0. 143R' ) * 

(T,p)  =  Pg(T,p)  +  vXT^p(AS’-0.143R')* 

where, 

AS'  =  entropy  of  melting/W^, 

V  =  (T-T  (o)}/{Tj^(p)  -  Tg(p)}, 

AT  =  Tj^(p)  -  Tg(p), 

\  =  {Tg(p)  +  T  (p)}/2, 

Tg(p)  =  solidus  temperature, 

Tjj^(p)  =  liquidus  temperature, 

X  =  coefficient  of  the  Lindemann  Law, 

=  d(ln  T  )/d(ln  p) 
tn 

In  reference  19,  the  AS  and  R  terras  appear  without  the  "primes,"  but 

the  primed  quantities  are  necessary  for  extrapolation  to  the 
solid  state,  and  because  of  the  stated  "scaled  entropy  of  melting." 


b.  Liquid-Phase  Region 


Ej^d.p)  =  EgCT.p)  +  Tg  AS'  +  0.143R'  AT/2 
+  lil  A(T)  , 

2a 


Pj^(T,p)  =  Pgd.p)  +  ATn,  pAS’  +  3R'-\  pTm  A(T) , 

2a 


where. 


A(T)  =  IniCST/Tm)  +  D  -  (l+a)ln(l+a) 

-  a{(T/Tn,)  -  1} 

a  *  fitting  parameter  for  liquid  specific  heat;  reference  19 
uses  a  =  0.1. 


c.  Vapor-Phase  Region 

E^(T,p)  =  Es(T,p)  +  3R’T/2  -  2R'T^p{(2-n)  /(1-n)^}  ^  -  ayp, 

dT 

P^(T,p)  =  R'Tp{(l+n+n2-n3)/(l-n.)2}-  ayp^, 

where, 

n  =  oVj,(T), 

=  vapor-exclusion  volume, 

Sy  =  coefficient  of  the  attractive  potential, 
dVb/dT  =  -Vb/(4T). 

Obviously,  the  above  equations  of  state  cannot  be  used  until  a  large 
amount  of  information  is  known  about  a  material.  Also  reference  19  presents 
a  program  for  estimating  the  critical  constants  associated  with  the  liquid- 
vapor  mixed-phase  regions,  and  these  constants  are  required  in  order  to  use 
the  previous  equations  for  a  general  mixed-phase  condition. 

We  now  need  to  develop  a  constitutive  equation  which  is  applicable  to 
both  solid  and  fluid  behavior. 


B.  Constitutive  Equation  for  "Hygrosteric”  Materials 


An  excellent  review  of  constitutive  equations  for  "plastic"  solid- 
behavior  is  given  in  reference  15.  However,  irrespective^®  of  the  temper¬ 
ature  dependency  of  the  material  parameters,  these  equations  do  not  take 
the  classical  Newton-Cauchy-Poisson  fluid-equation  form 


o  =  -pi  +  ,\(tr  A)1  +  2\if  A 


where , 


0  =  stress  tensor. 


I  =  identity  matrix. 


p  *  hydrostatic  pressure, 

A  =  matrix  with  components  A^^  =  jOV^/aXj  +  3V^/3X^), 

=  material  parameters. 

This  poses  a  problem  since  we  would  like  a  constitutive  equation 
which  continuously  maps  the  material  behavior  between  the  solid  and  liquid 
states.  Noll^^  discusses  the  concept  of  a  "hygrosteric”  material  which  can 
take  either  solid  or  fluid-like  behaviors.  Depending  upon  the  functional 
assumptions,  the  hygrosteric  constitutive  equation  for  an  isotropic 
material  can  be  written  as  a  general  function^^  of  the  strain  invariants 
dl,  I2  and  I3),  and  various  powers  of  (tr  A).  Any  general  function  of 
this  form  satisfies  all  the  restrictive  conditions^®"^^  associated  with 
the  mathematical  consistency  of  continuum  mechanics.  Also,  an  acceptable 
form  of  the  constitutive  equation  involves  the  decomposition  of  stress 

0  into  a  volumetric  component  pi,  and  the  so-called  "extra"  stresses  S. 

Of  all  the  equations^ ^  that  satisfy  the  above  general  forms,  the 
author  selected  the  simplest,  i.e., 

=  {-P(T,p)  +  K(T)(/r3-l)*  -  3K(T)a(T)(T-TD)}  6ij  +  s^j  (5) 


For  small  strains,  (/T3-l)-»-ei  i+e22+^33-^vi, 

C.  iT^uesdell,  "Elastiaity  and  Fluid  Dyncjnias,"  J.  Rational  Meahanics 
and  Analysis,  1,  126-262,  1952,  and  Vol  3,  693-611,  1953. 

■>1 

Walter  Noll,'  "On  the  Continuity  of  the  Solid  and  Fluid  States,  "  J. 
Rational  Mechanics  and  Analysis,  4,  3-81,  1955. 

)  o  ■  ■  -  —  - 

'“Walter  Noll,  "A  Mathematical  Theory  of  the  Mechanical  Behavior  of 

Continuous  Media,  "  Archives  Rational  Meah.  Analysis,  2,  198-226,  195c 
13  _ 

Proceedings  of  the  International  Conference  on  Constitutive  laDs  for 
Engineering  Materials :  Theory  and  .Application,  Tucson,  Arizona, 

Jan,  1983. 
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where  the  terms  refer  to  the  total  volumetric  component,  and  P(T,p) 

represents  the  EOS  pressure  component.  Equation  (5)  passes  in  the  limit 
(T-*-Tjj^(p) )  to  Equation  (4)  if  the  extra,  or  deviatoric,  stresses  s^j  depend 

upon  the  velocity  gradient  Aij  when  ’P+T^Cp).  This  can  be  accomplished  by 
using  a  modification  of  Perzyna's  constitutive  equations^**for  rate-sensitive 
materials  and  Zienkiewicz  and  co-workers  used  this  approach  in  references 
3  and  25. 

Perzyna  introduced  a  plasticity  function  ‘1>(F)  which  satisfies  the 
conditions 

(F)  =  0  when  F  <  0 
<j(F)  /  0  when  F>  0 

where  F  =  Prager's  dynamic  overstress-function 


=  ('yig  )/K(T)  -  1 

(2) 

I  =  second  invariant  of  the  extra  stress-tensor 
s 

=  ^  s  .  .  s .  .  , 
ij 

K(T)  ="yield  stress"in  simple  shear. 

Note  that  K(T)  can  also  be  a  function  of  hydrostatic  pressure,  effective 

(2) 

deviatoric  strain  Ig  =  '■-i  6^4  fiii,*  and/or  effective  deviatoric  strain- 

(2) 

rate  I.  =  4  e. .  e. . . 
e  "  ij  ij 

Cristescu  and  Predeleanu^^  manipulated  Perzyna's  equations  into  the 
following  form  for  loading  conditions  +  t  5<J>/5T}  >  0)  and  an 

"elastic  viscoplastic"  body. 


3y  definition,  e^ 4  =  e..  -  5..  (/la  -  1)11. 

ii  11 


?.  ?erzyna,  "The  Constitutive  Equation  for  Rate  Sensitive  Plastic 
y.azerials,"  Quarterly  Aopl.  Math.,  XX,  No.  4,  321-322,  1963. 

C.  Zienkievicz ,  ?.  C.  Jain  and  E.  Onate ,  "Plow  of  Solids  SuriKc 
Porring  and  Extrusion:  Some  Aspects  of  Numerical  Solutions,”  Inz. 
J.  Solids  Szructures ,  14,  15-33,  1973. 


(7) 


s  .  »  (1  +  6.  )*G  e.  .  if  J  <  K(T)^, 

ij 


+ 


\ff 


if  J  >  K(T)’, 


where  G  =  shear  modulus, 

J  =  is  {(1  +  V^^®ij  -  ^ 

eii^  =  s. ./(I  +  6. .)*G, 
iJ  ij  ij 

®ii^  =  values  of  e..  satisfying  J  =  K(T)^, 
J  ij 

Grj,  =  strain-hardening  shear  modulus 


and  viscosity  is  a  non-linear  function  of  stress,  strain,  strain- 

rate  and  material  properties.  For  unloading,  the  stresses  must  satisfy 

•  « 

the  inequality  ^ TS^/ST  <  0),  whereas  relaxation  implies  that 

^^ij^ij  T34>/3T  =  0). 

Equation  (7)  satisfies  both  the  solid  and  fluid  equations  if  the 
effective  viscosity  is  written  (assuming  a  rule-of -mixtures)  as, 


n 


eff 


f +  f  n  ^  , 

^  ^  ®  ®  K(T)<t(F) 


(8) 


where  f^ 
f 

s 

F 


volume  fraction  of  fluid, 
volume  fraction  of  solid, 

(/T /K)  -  1, 


ri-  *  fluid-viscosity,^®  which  can  be  a  function  of  temperature, 

(2) 

effective  deviation  strain  I  and  effective  deviatoric 

(2)  ® 

strain-rate  !•  , 

e 

n  =  solid-viscosity,^^  which  can  be  a  function  of  temperature 
and  the  strain-hardening  function  J, 


Because  of  the  strain  t..  definition,  (1  +  5..)  appears  instead  of  the 
usual  factor  of  2. 


'“'5.  L.  Sinkhoviah,  A.  S.  Romanov,  and  L.  I.  lonochkvna,  "Riidnoa'^nzmto 
Stability  of  a  Generalized  Nonlinear  Visco-Rlastia  Fluia  i'r.aer  Flano 
Gradient  Flow,”  translated  from  lyizhenemo-Fizioheskii  Zhuma  I ,  nc,  'S  . 
t'0-64,  2933,  Flenum  Publishing  Carp.,  1934, 


30 


limit  K(T)  -*•  0, 

T  ^  T^(p) 

limit  {fg/(K(T) $(?))}  0 

T  ^  Tj_(p) 

N^ce  that  the  above  second  limit  is  satisfied  for  any  of  the  non-linear 
?  (F-*F)  functions  proposed  by  Perzyna.  Also,  constitutive  equation  (7) 
provides  for  both  rate-independent  and  rate-dependent  material  behavior. 
Hence,  with  an  appropriate  selection  of  material  constants,-^  one  should 
be  able  to  approximate  both  instantaneous  and  non-instantaneous  wave- 
propagation  phenomena.^® 

Now  that  all  the  items  contained  within  the  energy  equation  have  been 
defined,  we  can  proceed  with  the  development  of  the  numerical  solution 
scheme . 


IV.  DEVELOPMENT  OF  THE  LINEAR  MOMENTUM  EQUATION  FOR  THE  NUMERICAL  SOLUTION 


As  previously  discussed,  we  want  to  formulate  the  basic  equations  with 
respect  to  the  FE  approximation.  This  approximation  involves  discretizing 
the  geometry  (spatial  variables)  to  obtain  a  set  of  time-dependent  equations 
and  then  solving  this  set  using  some  other  approximation. 


A.  Discretization  of  the  Spatial  Variables 


In  order  to  discretize  the  linear  momentum  equation  (B.4)  the  previous 
strain-displacement  and  constitutive  equations  must  be  introduced  for  their 
corresponding  terms.  Furthermore,  since  the  deviatoric  behavior  changes 
when  J>K(T)^,  two  different  sets  of  equations  must  be  developed;  one  for 
the  elastic  behavior  and  another  for  the  fluid-type  of  behavior. 


1.  Elastic  Behavior 


The  stress  vectors  [0^]  can  be  written  as  the  following  for  i  =  1,  2 
and  3  respectively. 


[cj  =  [D  J[e] 


(9) 


3hushan  j.nd  V.  E.  ■Jahsnayi,  "Propagation  of  Weak  Waves 


',n 


Plastic  and  Elast: 


Solids  vith  Interfaces , 


^nz 


Lozl  =  ■ 


[as]  =  [Dj[c]  - 


0 

P 

,0. 

0 

0 

LPJ 


where  p  =  P(T,p)  +  3K(T)a(T) (T-Tq) , 


[Dx]  = 


r(X+2G)  A 
0  0 


A  0  0  0  1 

0  G  0  0 

0  0  G  QJ 


[Dy]  =  ro 


0  0 

A  (A+2G)  A 
0  0 


LO 

0 

0 

A 


0 

0 


0 

0 


G  0  0l 
0  0  0 
0  0  G 

0  G  0 
0  0  G 


A  (A4-2G)  0  0  0 


A  *  vE/d+v)  (l-2v)  , 

V  =  Poisson's  ratio, 

E  =  Young's  modulus, 

[e]  =[£  e  e  c  e  e 

XX  yy  zz  xy  xz  yz-‘ 


(10) 


(11) 


Note  that  the  modulus  can  be  a  function  of  strain,  i.e.,  non-linear,  as  long 
as  there  is  no  hysteresis  upon  unloading.  Also,  the  bulk  behavior  Oe=K(T) 
-1)  has  been  included  into  the  first  matrix-product  of  equations  (9) 
through  (11) • 

Equation  (1)  can  be  written  in  matrix  form  as  the  following, 

[e]  =  [l][u]  +  [nLU]  /2  ( 12) 


* 

For  anisotropic  elastic  behavior,  these  matrices  can  be  rewritten  in 
terms  of  the  elastic  coefficients. 
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where  [l] 


[NLU] 


3/  3x 

0 

0 

0 

3/3y 

0 

0 

0 

3/3z 

3/3v 

3/3x 

0 

3/32 

0 

3/  3x 

0 

3  /  o  z 

3/3y 

([nlx. 

[u] 

([nlx. 

[U])^ 

([nly. 

[u] 

([nly. 

[U]) 

([NLZ. 

[U] 

([nlz] 

[U]) 

([nlx. 

[U] 

([nly; 

[U]) 

([nlx] 

[U] 

([nlz; 

[U]) 

>* 

[U] 

([nlz; 

[U]) 

[NLX]  =  (3/3x)  [I], 
[nly]  =  (3/3y)  [l], 
[NLZ]  =  (3/3z)  [l]. 


[l]  =  3x3  Identity  matrix. 


[u]  =  [u,v,w]T. 


The  last  remaining  step  is  to  relate  the  displacements  [u]  to  the 
material  velocities  [v]  using  the  kinematic  relationship, 


[U] 


4^V:(x,y 


uo(x,y,z) 


V2(x,y,2,T) 

to 

/'^{V3(x,y.z,T) 


+  g^Ct) }dT  +  VQ(x,y,z) 
+  g^(T)}dT  +  wg(x.y,z) 


where  =  li,  v,  w  respectively  for  i  =  1,  2,  3  ug,  vg  and  wg  represent 
the  initial  deformation  at  time  tg  and  g^,  g^  and  g^  represent  the  time- 
dependent  motion  associated  with  R(t)  and  u)(t)  . 


Using  the  ab^ve  equations  in  the  linear  momentum  equation  (B.-*),  and 
rearranging  all  the  non-linear  terms  to  the  right-hand  side  ot  the  cq'aation, 
yields  the  following  non-linear  integral-differential  equation, 

(3W  [D  ]  [D  ])[L][U]dVol 

dX  ^  dV  «  JZ  ^ 


L 
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=  -  /  pw  [FgldVol  -/  pW  [v]  dVol  +/  W  [n][ST]dS 
/v  /  V  ■'s 

+  fiiv  ^’)(p  ^)(P  dVol  +j\i  [v]p[Hl[v]dS 

V  3y  3z 

11 

-  S  [NL(K)]  , 

k:=i 

where  NL(1)  = /jW  [d  ][nT-U]  dVol, 

J  ^  ^ 

V 

NL(2)  =/iW  [d  ][NLU]  dVol, 

3y  y 

h’L(3)  =/jW  [d  ][NLU]  dVol, 

32  ^ 

V 

N’L(4)  =/jW  ra[D^][L][u]  dVol, 

N'L(5)  =/jW  Ca][d  ][l][u]  dVol, 

V  ^ 

NL(6)  =/j;  [a][d^][l][u]  dVoi, 

J  j2 

V 

NL(7)  =  -/[!][  (p  3W)(p  ^)(p  c^)]^dVol, 

y  9  X  9  y  9  2 

NL(8)  =  -yw[A][n][T]dS, 

NL(9)  =/3W  [a][D  ][NLU]  dVol, 

ML(IO)  =/9W  [a][d  ][NLU]  dVol, 

NL(ll)  =/iW  [a][D^][NLU]  dVol. 


Note  chat,  even  when  all  the  non-linear  terms  [nl]  are  identically  zero, 
equation  (14)  is  still  implicitly  non-linear  because  the  pressure  term  p 
depends  upon  the  density  p,  which  in  turn  depends  upon  the  deformations  [u] 
through  the  relationship  p  =  p^/yT^. 
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Based  upon  the  above  equation,  the  finite-element  spatial  discretization 
is  developed  in  Appendix  E-1.  We  can  now  consider  the  fluid-type  of  behavior. 


2 .  Fluid-Type  Behavior 

For  this  case,  previous  equations  (12)  and  (13)  are  still  valid,  and 
only  the  stress  vectors  must  be  redefined.  Furthermore,  this  definition 

depends  upon  whether  the  material  is  undergoing  loading,  unloading,  or 
relaxation.  For  unloading,  one  forces  the  stress  and  strain  to  follow  some 
function  of  the  elastic  slope  (usually  the  so-called  initial  tangent-modulus 
Eq)  up  to  some  predefined  set  of  stress  and  strains,  after  which  the  same 
constitutive  equations  are  used.  Therefore,  this  can  be  readily  handled  by 
the  logic  of  the  program  and  doesn’t  require  a  new  set  of  equations.  The 
relaxation  behavior  is  predicted  by  solving  a  specific  differential  equation, 
and  hence  requires  its  own  solution  scheme.  However,  this  scheme  will  not 
be  developed  for  the  current  project.  Therefore,  only  the  specific  set  of 
equations  associated  with  loading  will  be  developed  for  the  fluid-type  of 
behavior . 


The  corresponding 


stress  vectors  [oj^]  can  be  written  as. 


[o;]  =  [GX]([e]  -  [s^]) 

+  [ETAX]([e]  - 

[oj]  ■  [GY]([e]  -  [e'’]) 

tETAY]([4]  -  ti^]) 


Ccj]  -  [GZ]([e]  - 

+  [ETAz]([e]  -  [e^]) 


(P+Pe) 
0  ■ 

0 


(P+Pg) 

0 


0 

(P+Pe) 


where  [GX]  = 


[GY]  = 


(2Gt^)  0  0 

0  0  0 

.  0  0  0 


.0  0  0 
0  (2Gj)  0 

.  0  0  0 


0  0  o'], 

G.p  0  0 

0  G.j^  0  _ 

G.p  0  0, 

0  0  0 

0  0  G'j*  - 


(15) 


(16) 


(17) 
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[GZ]  = 


0  0  0  0  Gq 

0  0  0  0  0 

0  0  (2G,j.)  0  0 


[etax]  =r(2n^^^)  0  0  0  0 


[ETAY]  = 


0  n  ..  0 

eff 

0  On, 


0  0  n  ..  0 

eff 


0  (2n^^^)  0  0 


0 

0 

0 

0 

0 

0 

0 

'^eff 

0 

0 

0 

0 

0 

(2n^ff)  0 

0 

yy  ’ 

®zz  ’ 

e  , 
xy 

\z’ 

e 

yz 

E 

yy’ 

E 

e  , 
zz 

E 

^y’ 

E 

®xz  ’ 

£ 

e 

yz 

E 

yy’ 

.E 

®zz  ’ 

.E 

®xy’ 

•  E 

^xz  ’ 

.E 

e 

yz 

'eff  J 
0  1. 


r.E-i  r.E  .E  .E  .E  .E  .E  -iT 
LeJ=Le  ,e  ,e  ,e  ,e  ,e  J. 

^  xx’  yy’  zz  xy  xz  yz^ 

E  •  E 

Noce  that  the  individual  terms  e..  and  e..  are  determined  from  the 
previous  definitions. 

The  deviatoric  strain-rates  e..  are  related  to  the  strain-rates  t.. 


[e]  =  [e]  -  [0  0  0  f  f  f]^  (18) 

where  f  =  l/3(\/T3-l)  and  f  =  3f/3t.  The  term  VT3  is  the  determinant  of  the 
matrix  [a] ,  and 

'(l+3u/.3x)  0  0  1(0  )  Ou/cy)  (3u/3z) 

[a]  =  0  (l+3v/3y)  0  +  (3v/3x)(  0  )(3v/3z)  (19) 

0  0  (l+9w/3z)_  (3w/3x)  (3w/3y)  (  0  ) 


The  time  derivative  of  this  determinant  can  be  written  as  the  following 
determinant , 


3  C 


1  - 

=  3V 

3V 

av 

!  _ X 

X 

X 

3x 

3y 

3z 

' 

3V 

3V 

3V 

:  _ y 

_Z 

_ y 

!  3x 

‘ 

3y 

3z 

i  3V 

3V 

3V 

!  _£ 

z 

z 

i  3x 

3y 

3z 

(20) 


where  the  order  of  differentiation  was  interchanged  and  equation  (13)  was 
used  for  [u]. 


The  strain-rate  matrix  [e]  is  given  by  the  following  equation, 

[c]  =  Cl][v]  -h  [nLU] 


(21) 


where 

[nlu] 

6x 

2(^ 

V 

3  X  ^ 

3v 

V 
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V 
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,  3y 

i_n 

+ 

5w 

ay  ( 

jy 

3y  ; 

ay 

3y  3y 

“  32 

32 

32 

32 

32 

v.^ 

V 

V 

V 

V 

;,u 

an  a  X 

in  +  in  i-z 

+  1- 

2 

+  - 

''ax' 

ay 

3- 

s  sy 

ax 

■  ay  a 

X  ay 

3x 

iy 

V 

V 

V 

V 

V 

in 

.3  X 

in  + 

In 

i_ii  4 

.  3  y 

in 

+  ^ 

+ 

3  z 

4- 

3w 

32  ^ 

3x 

32 

3x 

32 

3x 

3z 

3x 

32 

3x 

V 

V 

V 

V 

in 

3w 

j.  3  X 

3v 

i_Ji  ^ 

.  3_n 

3v 

+  in 

i-X  + 

3  z 

+ 

^3y 

3z 

3y 

3z 

3y 

3z 

3y 

3z 

3y 

3z 

3y 

V 

",w  ;  z. 


3y 

.V 


Now,  using  equations  (15)  through  (21)  in  equation  (B.4)  and  manipula' 
ting  terms  yields  the  following  equation. 


f  sw 


^  [ETAX]  +  W  [ETAY]  +  W  [ETAZ])  [l][v]  dVol 
3x  3y  32 


=  -/ pW  [Fg]  dVol  -f  pW  [V]  dVol  +/ [n][ST] 
A,’.  A;  A 


dS 


+7[((P+Pe)£W'  ((P+Pe)3W)((p+Pe)^)]'  dVol 
■V  3x  3y  3z 

+  /"  W  [v]  0  [n][v]  dS.  -  /  (3W  [GX]  +  3W  [gy] 


(22) 


+  ^  CGZ])([l][u]  -  [h]  -  [e^])  dVol 
3z 


+  Ao^  [ETAX]  +  ^  [ETAY]  +  ^  [ETAZ])([e^]  +  [h]) 

J  3x  3y  3z 

-'v 


dVol 


=  -E  [nLF(K)], 

K=1 

where  NLF(l)  =  ^  [GX][nLU]  dVol , 

V 

NLF(2)  =,  ^  [GY][NLU]  dVol, 


V 


3y 


NLF(3)  =,  ^  [GZ][nLU]  dVol, 

NLF(4)  ^  [ETAX][NLU]  dVol, 

3  X 

V 

^^LF(5)  =,  3w  [eTAY][nLU]  dVol, 

NLF(6)  =  3W  [ETAZJCnLU]  dVol, 

3  2 

V 

NLF(7)  =.  IW  [a][GX]([l3[u]  -  [h]  -  [e^])  dVol, 

NLF(8)  =  ^  [a][gy]([l][u]  -  [h]  -  [e^])  dVol, 

_ 

NLF(9)  =,  iW  [a][GZ]([l][u]  -  [h]  -  [e^])  dVol , 

NLF(IO)  =  -/[A][((p+Pg)3W)((p+Pe)iW)((p+Pe)3W)]^  dVol, 


3x 


NLF(ll)  =  W  [n][ST]  dS, 

N'LFliZ)  =  .  3W  [a][GX][NLU]  dVol, 

'  3x 
.  V 

.VLF(13)  =  [a][GY][nLU]  dVol, 

'  5y 

V 

N'LF(14)  =  3W  [a][GZ][NLU]  dVol, 


3y 


3z 
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B.  Discrecization  of  the  Temporal  Variable 

The  equations  given  in  Appendix  E  represent  a  set  of  non-linear,  first 
order,  differential  equations  (actually  an  integral-differential  equation 
for  the  elastic  behavior)  for  the  unknown  material  velocities  [ v] .  Hence, 
an  additional  approximation  must  be  introduced  in  order  to  obtain  a  solution. 
Numerous  FE  programs  use  finite-differences  to  obtain  a  discretization  of 
the  temporal  variable,  with  the  explicit  forward-differences  being  the  most 
common  technique.  A  few  programs  use  the  implicit  Crank-Nicolson  (central 
difference)  scheme,  which  has  second-order  accuracy  with  respect  to  trunca¬ 
tion  errors  in  a  linear  problem.  Baker**  uses  the  Method  of  Weighted 
Residuals  (the  FE  method)  and  finds  the  approximate  solutions  to  be  as  good. 
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or  better  than,  the  solutions  from  the  Crank-Nicolson  method.  A  similar 
observation  was  made  in  reference  28  based  upon  simple  first  and  second 
order  equations.  Therefore,  the  author  decided  to  use  the  FE  method  for 
the  temporal  discretization. 

Depending  upon  the  choice  of  the  time-weighting  function,  a  wide 
variety  of  time-marching  equations  can  be  obtained,  with  some  equations 
yielding  better  approximations  than  others  (see  references  8  and  27).  A 
variety  of  final  equations  are  presented  in  Appendix  F,  and  during  the 
second  year  of  research  the  author  will  attempt  to  ascertain  which  equations 
yield  the  best  approximation. 


V.  DEVELOPMENT  OF  THE  ENERGY  EQUATION 
FOR  THE  NUMERICAL  SOLUTION 

The  energy  equation  (D.4)  must  be  rewritten  in  terms  of  the  unknowns 
and  T  by  using  the  EOS,  strain-displacement,  and  constitutive  equations. 
Furthermore,  the  equation  takes  different  forms  depending  upon  both  the 
thermodynamic  phase  and  the  deviatoric  stress  behavior,  i.e.,  whether 

J<K(T)^,  or  >K(T)^. 

Hence,  in  order  to  spatially  discretize  equation  (D.4),  several  new 
terms  must  be  explicitly  developed.  The  first  term  involves  the  time 
derivative  of  the  thermal  part  of  the  EOS.  This  part  can  be  written  as 
B(T)  +  B^(T,|/T^)  T  +  B2(T,/r2)  T^,  where  the  B  coefficients  depend  upon  the 

the  thermodynamic  phase  of  the  material.  Therefore,  the  time  derivative 
becomes , 

9  (thermal  part  of  EOS)  =  ^  ^  +  3B j  jT  T 
3t  dT  3t  3T  9t 


+  T  +  B^  ar 

d/Tj  dt 

+  iB^  ayT^  T^  +  2B2 
3t 

which  has  the  general  form 


+  38-  3T  T^ 
3T  3t 


T  3T 
3t’ 

of  B^(T)  +  B^CT)  T. 


(23) 


‘^C.O.  Zienkiewiaa  and  K.  Morgan, 
Chapter  ? ,  Cohn  Wiley  and  Sons, 
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Another  new  term  is  the  product  [e]  [o..]  of  strain-rate  and  stress. 

However,  the  strain-rate  is  discretized  in  Appendix  E  and  this  product  can 
be  readily  evaluated  by  using  a  slightly  different  form  of  equations  (9)  - 
(11)  and  equations  (15)  -  (17).  Therefore,  using  equation  (23)  and  the 
same  spatial  and  temporal  discretizations  as  used  in  the  linear  momentum 
equation,  the  energy  equation  can  be  fully  discretized  and  this  development 
is  presented  in  Appendix  G. 


VI  SUMMARY 


In  order  to  develop  a  predictive  2-D  model  for  penetration  problems, 
a  self-consistent  set  of  equations  was  developed.  This  set  of  equations 
includes  the  concept  of  polymorphic  and  thermodynamic  phase  changes,  and  uses 
a  corrected  BRL  version  of  the  GRAY  three-phase  equation-of-state .  Since  the 
critical  temperatures  and  pressures  have  already  been  established  by  the  BRL 
for  18  different  metals,  existing  data  and  FORTRAN  coding  can  be  incorporated 
into  the  new  model.  The  polymorphic  changes  are  incorporated  via  existing 
Hugoniot  data  which  exhibit  volumetric  discontinuity  at  various  pressures. 

In  addition  to  the  equations-of-state,  the  new  model  also  considers 
a  constitutive  equation  involving  both  the  scalar  pressure  (which  is  depend¬ 
ent  upon  volume  change  and  temperature)  and  a  deviatoric  stress.  This 
deviatoric  stress  contains  both  strain  and  strain-rate  contributions  with 
temperature  dependent  material  properties.  Furthermore,  this  deviatoric 
stress  maps  continuously  between  the  solid  and  fluid  (liquid  and  vapor) 
states  of  the  material.  Also,  the  deviatoric  stress  depends  upon  the  entire 
history  of  the  velocities,  and  hence  the  model  could  be  applied  to 
temperature-dependent  viscoelastic  materials  such  as  plastics. 

Because  of  the  thermodynamic  phase-changes ,  the  conservation  of  mass 
equation  is  written  in  a  global  fashion  to  keep  track  of  the  mass  contained 
with  the  solid-liquid-vapor  phases.  Furthermore,  mass  can  be  convected  out 
of  one  phase  into  another  and  the  corresponding  momentum  and  energy  fluxes 
are  considered  in  the  conservation  equations. 

Since  large  extensions  and  rotations  of  material  can  readily  occur 
during  a  penetration  process,  Che  new  model  uses  the  strain-displacement 
relationships  which  include  the  second-order  terms.  Also,  the  finite 
deformation  definition  is  used  to  represent  volumetric  changes  within  a 
material.  Lastly,  to  account  for  a  possible  reorientation  of  the  stresses 
during  a  time-step,  the  nonlinear  forms  of  the  linear  and  angular  momentum 
equations  are  used. 

Since  temperature  is  one  of  the  essential  unknowns,  the  energy  equa¬ 
tion  is  introduced  along  with  the  possibility  of  conductive,  convective  and 
radiation  heat  fluxes,  internal  heat  generation  due  to  latent  heat  and 
chemical  reactions,  internal  kinetic  and  stored  strain  energies,  work  of 
inertial  forces  and  applied  surface  tractions,  and  internal  irreversible 
work.  Note  that  generalized  inertial  forces  are  included  in  all  of  the 
conservation  equations  so  that  problems  involving  bodies  with  translational 
acceleration,  and  rotational  velocities  and  accelerations,  can  be  considered 
for  analysis. 

The  above  resulting  set  of  differential  and  integral-differential  equa¬ 
tions  are  extremely  nonlinear  and  are  coupled  together.  Even  when  applied 
to  the  idealized  problem  of  a  stretching,  one-dimensional  shaped-charge  jet, 
the  resulting  set  of  equations  cannot  be  directly  integrated.  Hence,  the 
author  introduced  an  approximation  by  discretizing  the  equations  in  both 
time  and  space.  This  discretization  involved  the  method  of  weighted  residu¬ 
als  using  piecewise  continuous  weighting  and  trial  functions.  This  converts 


the  original  nonlinear  differential  equations  into  a  set  of  nonlinear 
algebraic  equations  which  must  be  evaluated  at  each  time-step.  This 
incremental  solution  involves  algebraic  recursion  relationships  which  depend 
upon  the  assumed  weighting  and  trial  functions.  A  wide  variety  of  recursion 
relationships  were  derived  and  are  presented  in  the  appendices. 

The  beauty  of  this  newly  developed  set  of  equations  is  that  almost  all 
nonlinearities  appear  as  separate  entities  within  the  equations.  Hence,  the 
various  nonlinear  contributions  may  be  successively  "turned  on"  by  either 
the  user  or  some  sort  of  computer  logic.  Also,  the  user  can  specify  regions 
with  different  nonlinear  contributions  which  the  computer  logic  can  bring 
"in  and  out"  for  the  analysis  of  the  different  regions.  Furthermore,  one, 
two,  and  three  dimensional  contributions  have  been  individually  identified 
along  with  their  contribution  to  the  order  of  the  matrices.  Therefore,  the 
2-D  formulation  can  be  readily  updated  to  3-D  type  of  problems  once  the 
appropriate  functions  and  matrices  are  modified. 

The  last  aspect  of  the  new  model  is  the  "creation  and  annihilation" 
concept  where  the  original  geometric  configuration  is  subdivided  into  a  set 
of  super-elements,  each  of  which  are  assumed  to  be  quasi-independent.  Each 
super-element  can  have  its  own  degree  of  nonlinearity,  with  the  vast  majority 
of  the  super-elements  being  governed  by  the  fully  linearized  momentum  and 
energy  equations.  Furthermore,  each  super-element  has  its  own  set  of 
discretized  equations  which  must  be  solved  in  an  incremental  (and  possibly 
iteratively)  fashion  subject  to  its  boundary  conditions.  However,  many  of 
these  boundary  conditions  correspond  to  the  common  interface  between  various 
super-elements.  At  these  common  interfaces,  the  velocities,  tractions, 
temperatures  and  heat-flux  must  match  between  the  interconnected  super¬ 
elements.  Therefore,  these  interface  conditions  must  be  matched  via  an 
iterative  fashion  at  certain  time-steps  during  the  time  interval.  However, 
it  appears  that  this  iterative  C  and  A  process  will  yield  an  approximate 
solution  faster  than  current  FE  or  hydrodynamic  approaches. 
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LIST  OF  SYMBOLS 


C^(T)  = 

dVoi  = 

E  = 

E(T,p)  = 

e .  .  = 
ij 

* 

f  = 


coefficienc  of  specific  heat  at  constant  volume, 
elemental  volume. 

Young's  modulus 

internal  energy  per  unit  mass, 

ij^^  component  of  the  deviatoric  strain-tensor, 

inertial  body-force, 

l/3(^-l). 


f  = 
e 


13  = 

(2) 

I.  = 
e 

C2) 

I  = 
e 

(2) 

I  = 

S 

K  = 

n  = 
?(T,o)  = 
p„(.)  - 


volume  fraction  of  reversible  "plastic"  work, 

volume  fraction  of  fluid  material, 

volume  fraction  of  irreversible  "plastic"  work, 

volume  fraction  of  solid  material, 

shear  modulus, 

strain-hardening  shear  modulus, 
third  strain-invariant  of  the  e  strain, 


4  e .  .  e  ,  ,  , 


'-5  e  .  .  e  .  .  , 
iJ  iJ 


h  s..  s. . , 

yield  stress  in  simple  shear, 

thermal  conductivity  in  the  principal  material  direction, 

outward  unit-vector  normal  to  a  surface, 

F.quation-of -state  (EOS)  pressure, 

Kugoniot  pressure. 


p  =  P(T,p)  +  3<(T)ci(T)(T-Tp), 
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heat  source  or  sink  per  unit  mass, 
specified  heat-flux  across  a  surface, 

global  position  vector  with  respect  to  inertial  coordinates. 
Universal  Gas  Constant/ Atomic  weight. 

Local  position  vector  with  respect  to  an  attached  rectangular 
cartesian  coordinate  system, 

surface  separating  two  different  material  phases  at  time  t, 
ij^^  component  of  the  "extra"  or  deviator ic  stress-tensor, 
temperature  K, 

Debye  temperature, 
liquidus  temperature, 
solidus  temperature, 
time , 

initial  time, 
shock  velocity, 
particle  velocity, 

volume  of  the  sub-domain  at  the  start  of  a  time-otep, 
material  velocity  in  the  i^*^  direction. 

Spatial  Weighting  function. 

Atomic  weight, 

temporal  weighting  function, 

unidirectional  coefficient  of  thermal  expansion, 

T^(p)  -  Tg(p), 
time-step  size, 

ij^^  component  of  the  strain  tensor. 


effective  material-viscosity,. 


elastic  bulk-modulus  coefficient. 

Lame's  constant, 
mass  density, 

ij^^  component  of  the  stress  tensor, 
time-variable, 

dimensionless  time-variable, 

Perzyna's  plasticity  function, 

rotational  velocity  about  the  i^^  axis  of  a  coordinate  system 
angular  acceleration  about  the  i  axis, 
rotational  angle, 

a  matrix  with  n-rows  and  m-columns. 


Integrating  over  the  volume  V  of  any  particular  phase,  and  considering 

possible  phase-boundaries  S  ,  the  mass  conservation  equation  can  be  written 

P 

as , 


y'p(t^)dVol  -y'p(t)dVol  = 


v(.„) 


V(t) 


V’nds  )dT 
(T) 


(A.l) 


where  p  =  mass  density, 

dVol  =  elemental  volume, 

n  =  outward  unit  vector  normal  to  the  phase-boundary, 
V  =  material  velocity  across  the  phase-boundary. 


Note  that  the  term  inside  the  parentheses  on  the  right  hand  side  of  the 
equation  (A.l)  is  the  instantaneous  mass-flux  across  the  phase-boundary  and 

will  be  denoted  as  M  (t) .  The  total  mass  of  the  system  is  conserved  and 

P 


can  be  written  in  the  global  sense  as 


M 


=  p^dVol  +  y  p^dVol  +  j'p^ 


dVol 


V3(t) 


V^(t) 


V  (t) 

V 


(A. 2) 


where  V  (t) 
s 

V,(t) 

V^(t) 


p.. 


volume  of  solid  at  time  t, 
volume  of  liquid  at  time  t, 
volume  of  vapor  at  time  t, 
solid  density,  and 
liquid  density. 


vapor  density. 


APPEND ly.  B 

CONSERVATION  OF  LINEAR  MOMENTLIl  EQUATIONS 


Based  upon  the  assumptions  discussed  in  the  text,  Novozhilov's  linear 
momentum  equations  can  be  written  in  the  following  elemental  form. 


(3.1) 


^  ■  ([A^]  [Oj^]i  +  [a^]  [a^lj  +  [a^]  [o^Jk) 

-  pFg^  =  p  +  (V^So)*  i=l,2,3 


where  V  =  ^  i  +  ^  j  +  3_  k, 
9x  3y  9z 


[A^]  =  row  matrix  aasnriated  with  the  ith  row  of  [a], 

[a]  =  (l+3u/9x)  (9u/9y)  (9u/9z) 

(9v/3x)  (l+9v/9y)  (9v/9z) 

(9w/3x)  (9w/9y)  (l+9w/9z) 

[a.]  =  column  matrix  associated  with  ith 

column  of  the  stress  matrix  [a], 

[a]  *  a  c  0  , 

XX  yx  zx 

a  0  a 

xy  yy  zy 

a  a  a 

xz  yz  zz 

u,  V,  w,  =  components  of  the  displacement  vector, 

R.  =  time  derivative  of  the  ith  component  of  the  global 

position-vector  R, 

=  ith  component  of  the  inertial  body-forces  per  unit  mass, 
a..  =  ijth  component  of  stress, 

=  u,  V  or  w  for  i  =  1,  2,  3  respectively. 

Applying  the  method  of  weighted  residuals,  equation  (B.l)  is  multi¬ 
plied  by  some  weighting  function  W (x,y,z) ,  and  the  resulting  equation  is 
integrated  over  the  entire  body,  i.e.. 


Jw  V-([a.]  [o,]i  +  [a.]  [a.]j  +  [A.]  raJk)dVol-  /  pWF„  .  dVol 

V(t)  1  1  1  2  1  3  Bi 

/pW  SV.dVol  -  f WV. (pV.n.)dS  , 

» T  /  .  N  O  t.  _  » 


i'B.2) 


'.n.)dS  , 


Sp(t) 


Term  added  by  the  author  to  represent  possible  mass-change  within  the 
elemental  volume. 


where  is  the  direction  cosine  for  the  ith  spatial  direction.  The  last 
term  represents  the  l,inear  momentum  flux  across  the  phase  boundary,  and 
corresponds  to  the  (V.9p/dt)  term  in  equation  (B.l)  coupled  with  the  con- 
servation-of-mass  equation  in  equation  (A.l). 


The  first  term  in  equation  (B.2)  has  the  form  (f  div  u)  which,  from 
calculus,  equals  (div  (fu)  -  u-gradf),  and  the  volume  integral  over  div(fu) 
can  be  replaced  by  a  surface  integral  using  the  Gauss  divergence  theorem. 
Using  these  facts  yields  the  following  weak-form  of  the  weighted  residuals. 


-  [o^]i  +  [A^] 

V(t) 


[02^j  +  [a2]k)-VWdVol 


•-  /w(Ca^]  Ca^]i  +  [A^]  [o2]j  +  [A^]  [ 

S(t) 


a^lk) -ndS 


-J'oV'F^AVol  =  /pW  aVi 
V(t)  'V(t)^‘^ 


,•  dVol 


yW  V.  (pV.n.)dS 
1  1  1 


S  (t) 
P 


(B.3) 


where  S  represents  the  surfaces  where  surface  tractions  t^^  are  specified. 

The  above  three  equations  can  be  written  as  a  single  matrix  equation  of 
the  following  form. 


3W  [a. j  +  3W  [o 
3x  3y 


1 

2^ 


V(t) 


+ 


3W  [a  ])dVol 
92 


V(t) 


dVol 


9  [7]dVol 


+  /pW  _^! 

J 

V(t)  ; 

+  /’(3W  [J][p  ]  +  9W 
y  9x  9y 


^«(Cn]  [ST]  +  [5]  [n]  [ST])dS 

(t) 

[3:][pJ  +  w  [S]  [a  :)dVol 
9z 


(B.4) 


P  [V] 
(t) 


[n] 


[V]dS 


60 


[n]  *  [n  n  n  ] 
X  y  2 


where  [n] 


[ST] 


r__-^ 

L^J 


Bi 


Q. 

1 


n  00 

X 


0  n  0 

* 

V 

0  0  n 

z 

(t  + 

t 

+ 

t 

XX 

yx 

zx 

(t  + 

t 

+ 

t 

yx 

yy 

yz 

(t  + 

t 

+ 

t 

zx 

yz 

zz 

(3u/9x) (3u/3y) (3u/32) 
(3v/3x) (3v/3y) (3v/3z) 
(3w/3x) (3w/3y) (3w/3z)  , 


coordinate  system  about 
ith  axis  of  the  local  fixed  coordinate  system, 


=  (0.  r,  -  0,  r.)  +  2(0.  V,  -  0,  V.)  + 

j  k  k  j'  ^  j  k  k  j 

(P..  r.  +  0,  r,  )  0.  -  (0^  +  P.^)  r., 

1  J  k  k  1  "  j  k  1 

=  rotational  velocity  of  local  material 


=  angular  acceleration  about  the  ith  axis, 


r^  =  ith  component  (r^,  r^  or  r^)  of  the  local  position  vector  r. 


Note  that  the  [a]  term  corresponds  to  finite  rotation  of  the  infini¬ 
tesimal  volume,  and  is  usually  assumed  as  zero.  The  first  term  in  the  body- 


force  Fg  expression  is  usually  called  the  "linear-acceleration"  body-force, 
the  second  term  is  the  "coriolis"  body-force,  the  third  term  is  the  classical 


The  elemental  form  of  the  angular  momentum  equation  is  developed  by 
taking  the  cross  product  between  the  position  vector  R  and  the  stress 
vectors, and  equating  this  to  the  time  derivative  of  (i?x(rav)).  Multiplying 
these  products  by  the  weighting  function,  using  the  fact  that  ^xv=0,  .'?=R+r, 
and,  R  is  independent  of  dVol,  integration  over  the  volume  yields  three 
equations  of  the  following  form. 


yw^r^V- ([Aj^]  [a^]i+[Aj^]  [o^jj+CA^]  Co^Jk) 

Ht) 

[a^]i+[A^]  Co2]j+CA^]  [a^]k)  |  dVol 

/p«(t.F3^-r^rjj)eVol  .  /'.■p<r^^-t^3V.)<iVol 

r 3t 


(C.l) 


minus  the  Weighted  Angular  Momentum  Flux  across  phase  boundaries,*  plus  one 
vector  equation  of  the  form  Rx  left-hand-side  minus  right-hand-side  of 
equation  (B.4).  Therefore,  when  the  linear  momentum  equation  (B.4)  is 
satisfied,  this  vector  equation  is  identically  equal  to  zero. 

Equation  (C.l)  can  be  rewritten  using  the  calculus  identity  involving 
fdivu  =  div(fu)-u*gradf ,  with  f=Wr.,  and  then  using  the  Gauss  divergence 
theorem  to  yield, 


j-  /([Aj^]  Cc^]i+[.A^]  [a2]j+[Aj^]  [a^]k)  •7(r^W)dVo] 
V(t) 

-^yrjW([Aj^]  [a^]i+[Aj^]  Cc^lj+CA^]  [o^jk) -ndS 
S(t) 

-  /  r  .pWF„,  dVol  -  /r.WpSV,  dVoll 


(C.2) 


-  ([Aj]  [o^]i+[A^]  Ca2]j+[Ajl  [a^]k)  •7(rj^W)dVol 


Coi]i+[Aj]  [a^^j+CA^]  f 

S(t) 

-  /r,  pWF„  .dVol  -  /  r,  Wp3V.  dVol . 


[o^ jk) -ndS 


+ith  Component  of  Weighted  Angular 
Momentum  flux  =  0 


*Note  that  this  quantity  corresponds  to  the  weighted  volume-integral  of  the 
j rx i  R+\r+.'.xr)Mp  1  terms. 


Comparison  of  the  terms  within  each  bracket  of  equation  (C.2)  with  the 
linear  momentum  equations  (B.3)  shows  that  each  terra  corresponds  to  a  linear 
momentum  equation  whose  integrand  is  multiplied  by  a  component  of  the 
position  vector.  Obviously,  if  the  linear  momentum  equation  is  satisfied  at 
all  points  within  the  body,  then  equations  of  the  form  of  equation  (C.2) 
are  also  satisfied.  However,  the  linear  momentum  equations  are  only  satis¬ 
fied  in  an  integral  sense,  and  therefore  the  angular  momentum  equations  will 
be  automatically  satisfied  only  for  certain  special  situations: 


a.  the  linear  momentum  integrand  is  independent  of  the  position 
component  r^ 

b.  the  linear  momentum  integrand  is  symmetrical  with  respect  to  the 
position  component  r.. 


The  conservation  of  energy  is  usually  stated  as  the  time-rate  of  change 
of  the  internal  energy  being  equal  to  the  sum  of  the  energy  flux  into  the 
system  plus  the  time-rate  of  change  of  the  heat  added  to  the  system  and  the 
work  done  on  the  system.  Based  upon  the  assumptions  mentioned  in  the  text, 
the  components  of  the  energy  equation  for  an  elemental  mass  can  be  written 
as 


_^(heat  added) 
3t 


_^(work  done) 
3t 


T 

=  pQdVol  -  (k  3T  n  +  k  3T  n  +  k  3T  n  )dS+f  [e]  [o..]dVol 

^  ^^37  "^31  P 

(D.  1) 

=  [Slf  [n]"^  ([II  +  [A]’^)|[R]+[V>[nxr]|dS  (D.2) 


_^(internal  energy)  = 
3t 


3  (thermal  energy)+f^[£  ]  ^  ^®ij^ 
3 1 


(D.3) 


+  9  (  P  ([R]'^+[V]'^+[f4xr]'^)  ([R]+[v]+[nxr])  >  dVol  . 

3t  2  ) 


In  equation  (D.l),  the  Q  term  is  positive  for  an  exothermic  heat  source 
per  unit  mass,  Fourier's  Law  is  assumed  to  hold  for  conductive  heat  flow, 


T  represents  the  absolute  temperature,  is  the  time-rate  of  change  of  the 


ijth  component  of  the  strain,  is  the  ijth  component  of  the  stress,  and 


f  is  the  volume  fraction  (0^  <.!)  for  non-elastic  deformations.  In 
P  ^  P 

equation  (D.2),  fl  is  the  rotational  velocity-vector  (Q^^+fi^+Q^)  ,  and  in 


equation  (D.3),  f  is  the  volume  fraction  (f  +f  =1)  for  elastic  deformations 
e  e  p 


Substituting  the  above  definitions  into  the  conservation  of  energy 
equation,  multiplying  by  a  weighting  function  W,  and  integrating  over  the 
entire  spatial  domain,  yields  the  following  equation, 


/ 


W  ;  p  3  (thermal  part  of  EOS)  +  f  [e]  [o  .1 
—  e'-  ' 


dVol 


3t 


V(c) 


^  •  -.T  -iT 

+lR]  *  weighted  linear  momentum  Eq  +  [o]  *  weighted  angular  momentum  Eq 


1 


3W[a,]  +  3Wra„]  +  3w[0-,]  +  pW([F„]  +  [v])  +  3WrA][0,  ]+  3U’[a][c, 


+|W[X'[a2]|dVol  - yivl'^W  jCn][ST]  +CS’][n]CST] j dS 

^  S(C) 

/•'■(pQ+f  [e]'^  +  /W(k  3T  n  +  k  3T  n  +  k  3T  n 

P  J  .XX-  X  yy-  y  zz-  . 


W(inCernal  energy) (pV*n)dS. 


Note  chat  the  second  and  third  terms  in  equation  (D.4)  are  identically 
zero  when  the  weighted  momentum  equations  are  satisfied.  The  integrands  of 
the  second  and  third  integral^  in  equation  (D.4)  involve  the  linear  momentum 
equation  post-multiplied  by  [V]'^,  and  the  last  integral  represents  the 
internal  energy  flux  and  includes  all  the  terms  involving  M  .  The  surface 
S  refers  Co  Che  surfaces  with  specified  heat-fluxes.  ^ 


In  order  to  apply  Che  C  and  A  process  to  the  linear  momentum,  we  must 
discretize  the  region  into  a  sub-set  of  NR  regions  or  superelements .  The 
unknown  [V]  within  each  of  these  regions  is  approximated  by  a  piece-wise 
continuous  function  of  the  form. 


[V]j  0  0  N,  0  0 

0  N,  0  0  0 


0  0  N,  0  0 


(E.l) 


where  NP(j)  =  number  of  nodal  points  in  the  jth  region, 

=  the  unknown  value  of  the  velocity  at  Che  ith  node  and 
in  Che  kth  direction, 

=  an  interpolation  function  in  terms  of  the  spatial 
variables . 

Hence,  all  Che  integrals  in  equations  (14)  and  (22)  become  a  sum  of 
integrals  over  each  region.  Note  that  the  spatial  variation  of  the  velocity 
is  now  expliciCy  given  by  Che  interpolation  function. 

Next,  the  weighting  function  is  assumed  to  be  approximated  by  a  similar 
piece-wise  continuous  function,  with  a  polynomial  N.(x,y,2)  associated  with 
each  node.  Since  equation  (14)  must  be  satisfied  for  each  discrete 
weighting  function  one  obtains  a  set  of  NP(j)*NDOF  simultaneous 
equations  for  each  region,  where  NDOF  refers  to  the  spatial  degrees-of- 
freedom.  Note  that  in  a  full  F.E.  approximation,  the  solution  for  all  the 

unknown  nodal  values  (DOF*^^P(j))  must  be  solved  simultaneously . 


1.  Elastic  Behavior 


The  [Lj[UJ  displacement  matrix  in  equation  (14)  can  be  written  as. 


[l][U]  =  [l1  [n] 

NFGXNDOF  OTOFxHT 


t-fA  t_ 

U  +  HV.  +  f  V.(T)d(T) 

L  °  ^  ^  JmTxI 


whe  re  HV . 

1 


history  of  the  velocity  at  the  ith  node  up  to  time  t. 


t 

/ 


V^(T)d(T) , 


to 


c 


o 

MT 


initial  nodal  displacements  at  time  t^, 


at 

NFG 


NOOF'  *  NP(j) , 

size  of  the  current  time-step. 


digit  which  depends  upon  Che  geometry  of  the  problem  being 
solved . 


Note  that  the  only  unknowns  in  equation  (E.2)  are  the  nodal  velocities  at 
Che  current  time-step  (the  initial  nodal  velocities  must  be  specified  at 


time  t^) 


The  two  remaining  velocity  terms  are  the  [NLU]  and  [A]  matrices  which 
can  be  written  as. 


FNLUI 


L 


[IV]^ 

[cxj" 

[cx] 

[IV] 

[CYl" 

[CY] 

[IV] 

[IV]" 

[cz]" 

[cz] 

[IV] 

[IV]" 

[CX]" 

[CY] 

[IV] 

[IV]" 

[cx]" 

[CZ] 

[IV] 

[IV]" 

[CY]" 

[CZ] 

[IV] 

(E.3) 


6x1 


and 


0  (UYIV)  (UZIV) 

(VXIV)  0  (VZIV) 
(VDCIV)  (WYIV)  0 

t+Ac 


(E.4) 


where  [IV] 


Fu  +HV.  +  /  V.  (T)dT] 

^  ^  -I 

-  =  F  — 1  °  °  —  0  0...  0  d 

]=FoiNi  003^0...0  3N^p  0 

L  3X  3x  3>:‘  J 


IT 


In  equation  (E.l),  NDOF  was  set  to  a  value  of  3. 


[ra]  =  [00^2  00^2  •••  Q09%p] 
3x  3x  3x 


[UZ]  =  [3^  00  3N2  00  ...  ^^^.pOO] 

3z  3z  oz 

[VZ]  =  [0  3N,  00  3N,  0  ...  0  3N  01 

■ 

[wz]  =  [00  3Ni  00  3N„  ...  00  3N  1 

[cx]  = 


[CY]  = 


[CZ]  = 


UYIV  = 

and 

WYIV  =  [WY]  [IV].  Note  that  NDOF  was  set  equal  to  its  maximum  value  of 
3  in  the  above  equations. 

The  remaining  terms  in  equation  (14)  that  must  be  spatially  discretized 
are  those  involving  the  weighting  function  W.  Using  the  previously  defined 
interpolation  functions,  the  terms  involving  the  weighting  function  become 
the  following, 


3  W  r  D  1  ~  ^  cx  1  ^  ^ 

^  -  X'  "  -'mTxNDOF  ^  x-NDCFxNFG’ 

3x 

(E.5) 

3W  [D^]  =  [CYf  [D  ], 

( E  .  6 ) 

[D^]  =  [UZ]'^  [D  ], 

(E.7) 

dZ 

WlFg,;  =  lNj  l-Fg 

^E.S) 

oz 


3z 


J  3xMT 


[UX] 

[VX] 

[wx] 

[UY] 

[VY] 

[U-Y] 

[UZ] 

[VZ] 

'  3xMT 
[UY]  [IV], 


3jeiT 


^E.S) 


w  [V]  =  [wf  [N]  [V], 


(E.9) 


W  [n][ST]  =  [N]  [n][STJ 
W  [v]  =  [N]^  [N][V1 


(E.IO) 

(E.ll) 


where  [N]  =  Nj^  0  0  0  0 

0  Nj^  0  0  N,  0 

0  0  N,  0  0 


0  0 
0  0  N., 


Equations  (E.l)  through  (E.ll)  are  substituted  into  equation  (14)  to 
give  the  following  set  of  spatially  discretized  equations, 

^([CXj^  [D^]  +  [CY]^  [D  ]  +  [CZ]^  [d^]Kb]  [lv]dVol 
V(j) 

=-rp[N]'^[F„]dVol  -  rpCNj"^  [N][v]dVol  + 

V(j)  V(j) 

J*  [N]''^[n][ST]ds  +  /p([UX]''^+iVY]''^+[WZ]'’^)  ^^3  dVol  (E.12) 

S(j)  V(j) 

r  T  -  -  11 

+  /[N]  [N][v]p[n][N][V]ds  -I  [NL  (k) ] , 

Sp(j) 


where, 


[B]  -  ^^^NDOFkMT  ’ 

NL(1)  =  /  [Cxf  [D  ]  [NLU]  dVol, 

V(j) 

NL(2)  =  J  [CY]^  [D  ]  [NLU]  dVol, 

V(j)  ^ 

NL(3)  =  .  /  iCZ'7  [D  ]  [NLU]  dVol, 


NL(4)  =  /  [CX]  [a]  [D^]  [B]  [IV]  dVol, 

V(j)  ^ 

NL{5)  =  /  [CY]''^  [A]  [D  ]  [B]  [iV]  dVol, 


NL(6)  =  J  [CZ]  [a]  [Dy]  [B]  [IV]  dVol, 

V(j) 

NL(7)  =  -/  [a]  p  ([UX]^  +  [VT]'^  +  [WZ]^) 

V(j) 

NL(8)  =  -/  [n]"^  [a]  [n]  [ST]  ds, 

S(j) 

NL(9)  =  /CCX]'^  [a]  [D^]  [NLU]dVol, 

V(j) 

NL(IO)  =  /[CY]^  [A]  [Dy]  [NLU]  dVol, 

V(j) 

NL(ll)  =  fiCZ]^  [A]  [D^]  [NLU]  dVol, 

V(j) 

V(j)  =  Volume  of  the  jth  region, 

and 

S(j)  =  Surface  area  of  the  jth  region  with  specified  surface 
tractions- 

Once  the  volume  and  surface  integrations  are  performed,  Equation  (E.12) 
represents  a  set  of  non-linear,  first-order ,_differential  equations  with 
respect  to  the  MT  unknown  nodal  velocities  [v].  Usually,  the  integrations 
are  done  using  the  numerical  method  of  Gaussian  Quadrature,  with  the  number 
of  integration  points  depending  upon  the  interpolation  function  and  the  type 
of  problem- 

2.  Fluid-Type  Behavior 

In  addition  to  the  discretized  matrices  given  in  the  previous  section, 
equation  (22)  involves  several  other  velocity-dependent  matrices.  The  most 
complex  is  the  strain-rate  of  equation  (21)  which  involves  a  matrix  of 
products,  with  each  product  discretized  to  a  series  of  matrix  multiplica¬ 
tions,  i.e.  , 


NLU 

=  [(U^/WX) 

(UVWY) 

(UVWZ)  (U^WNY) 

(UVIvTCZ) 

(UVWYZ) ] 

UVWX 

=  2 1 (UXIV) 

(UXV) 

+ 

(VXIV) 

(VXV) 

+ 

(WXIV) 

(^«V)|  , 

UVWY 

=  2|(UYIV) 

(UYV) 

+ 

(VYIV) 

(VYV) 

+ 

(WYIV) 

(WYV)|  , 

UVWZ 

=  2|(UZIV) 

(UZV) 

+ 

(VZIV) 

(VZV) 

+ 

(WZIV) 

(WZV)|  , 

UVWXY  =  (UXIV)  (UYV)  +  (UXV)  (UYIV)  +  (VXIV)  (VYV)  +  (VXV)  (l^’IV)  + 
(U-XIV)  (lA’V)  +  (WXV)  (WYIV)  , 


UVWXZ  =  (UXIV)  (UZV)  +  (UXV)  (UZIV)  +  (VXIV)  (VZV)  +  (VXV)  (VZIV)  + 
(WXIV)  (WZV)  +  (WXV)  (WZIV), 


UVWYZ  =  (UYIV)  (UZV)  +  (UYV)  (UZIV)  +  (VYIV)  (VZV)  +  (VYV)  (VZIV)  + 
(WYIV)  (WZV)  +  (WYV)  (WZIV), 

UXV  =  [ux]  [V], 

and  _ 

WZV  =  [wz]  [V]. 

The  remaining  matrices  are  [h]  and  [H]  which  involve  the  velocity 
implicitly  through  the  VT^  term.  Using  equations  (19)  and  (E.4),  and 
previous  matrix  definitions,  the  following  determinant  is  obtained. 


(1  +  UXIV)  (UYIV)  (UZIV) 
(VXIV)  (1  +  VYIV)  (VZIV) 
(WXIV)  (WYIV)  (1  +  WZIV) 


(E.14) 


The  time-derivative  of  this  determinant  can  be  written  using  equation 
(20)  as, 


JyJu  =  (UXV)  (UYV)  (UZV) 
at  ^  (VXV)  (VYV)  (VZV) 
(WXV)  (WYV)  (WZV) 


(E.15) 


Using  the  above  discretized  matrices,  equation  (22)  can  be  written  as 
the  follow  set  of  equations, 

y  ([CX]’’^  [ETAX]  +  [CY]’’^  [ETAY]  +  [CZ]'^  [ETAZ])  [B][V]dVol 


V(j) 

=  -/  p 

[Nf 

rFg]dVol 

-/p[N] 

V(j) 

V(j) 

/[Nf 

[n] 

CST]ds  + 

/p+Pg) 

S(j) 

V(j) 

+  J [N]"  [N]  [V]  p  [n]  [N]  [V]dS 


Sp(t) 


-  [GX]  +  [CY]'^  [GY]  +  [CZ]’’-  [GZ])  ([B]  [IV] 

V(j) 


[H]  -  [e^])dVol 


+  y([CX]’^  [ETAX]  +  [CY]"^  [ETAY]  +  [CZ]’’’  [ETAZ]) 
V(j) 


(E.16) 


-Z  [NLF(k)]. 
k=l 

Note  that  in  equation  (E.16),  the  6x1  matrices  [e  ],  [e  ],  [H]  and  [h] 
involve  array  elements  which  are  the  matrix  products  of  all  the  correspond¬ 
ing  MX  nodal  values.  The  twenty-three  NLF  terms  are  not  rewritten  here 
since  the  forms  are  the  same  as  before  except  with  V->  V(j),  S-»-  S(j),  and 
equations  (E.5)  through  (E.ll)  being  substituted  for  the  corresponding 
terms • 
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In  both  equation  (E.12)  and  equation  (E.16),  the  spatial  variables 
have  been  integrated  out,  and  the  remaining  unknowns  are  the  MT  nodal 
velocities  as  some  continuous  function  of  time.  The  next  approximation 
is  to  replace  the  continuous  time-behavior  with  a  finite  set  of  times  for 
determining  the  nodal  velocities. 

1.  Elastic  Behavior 

Equation  (E.12)  can  be  written  in  the  standard  form  as 

a.  _ 

([KX]+rKY]+[KZ])[rv]=  -[FB]+[FT]  -[m][V]+[P]+[MF][v]-Z  [NL(k)],  (F.l 

k=l 


where , 


[KXl=  /[CX]  [D  ][B]dVol, 
V(j)  "" 

[KY]=  /[CY]’^[D  ][B]dVol, 


[KZ]=  /[CZ]  [D  ][B]dVol, 

V(j) 

[FB]=;p  [xfCF  ]dVol, 

V(j)  ® 

[FT]=  .^LN]'’’[N][ST]dS, 


Ch]=/0  [N]  [N]dVol, 

V(j) 

[P]=/P  ([UX]'^+[VY]^+[WZ]'^)  dVol, 

V(J) 

and,  [MF]=/  [N]'^([N][v][n])  [N]dS. 

Sp(j) 


We  now  assume  a  piecewise  non-zero,  time-dependent,  weighting  function 
WT(t),  and  approximate  the  time— dependent  nodal  velocities,  using  the 
following  two  equations,  respectively. 


WT(t)=  WT  (t), 
n 

n*0  . 


(F.2) 


and  [V(t)]  =ENT  (t)[V(n)], 
MTxl  " 


(F.3) 


where  is  the  total  number  of  time-steps  selected  for  a  particular  problem, 
[V(n)]  represents  the  MT  nodal  velocities  at  the  time  t  ,and  NT  (t)  is 


:;'.o  ir.cerpol:!!;  ion  fonciion  ror  c’nc  nch  t  ir.o-o  ler.on  C .  ^o:e  Co.at  Cr.a  ...^(t;) 

functions  are  non-cero  only  over  che  ntii  cir.e-eler’.onc .  Substituting 
equations  (r.2)  and  (r.3)  into  equation  (r.l)  yields  a  set  of  equations 
of  the  folLo'.ving  for:?. 


J'(_  nX  , In':' ^ KZ  j )  _  I (  t)  d  t  - 


/  {-[?'3]+tE'T]-rL?j  -  Z  ^\’L(k)  ]t'.-.T^(t)dt 
r  ’k.=  1 


(F.i) 


■(jLM:-T  (t)';NT  (t)  ]dc)[V(n)] 
"In  n 


(t)LNT  (c)]dc)[v(n)], 

n  n 


uhere  t  =  tir.e-snan  over  vhich  che  nth  voighting  function  is  non-zero, 
s 


A  wide  variety  of  equations  can  be  obtained  from  equation  (F.4) 
depending  upon  the  .assumptions  involving  VT  ,  N’T  and  t_.  Also,  the  terms 
inside  the  spatial  matrices  may  change  with'  time' due  to  either  changing 
boundary  conditions,  or  temperaruce-dependent  material  properties.  Since 
the  characteristic  thermal  time"  can  be  several  orders  of  magnitude  less 
chan  che  characteristic  mechanical  ti.me“^  tor  impact  problems,  temperature- 


dependent  material  properties  can  be  assumed  constant  during  time  steps 
small  in  comparison  to  the  characteristic  thermal  time.  The  boundary 
conditions  are  usually  matc'ned  by  using  sufficiently  small  time-steps 
associated  with  e."plicit  solution  schemes.  However,  for  implicit  schemes, 
which  use  larger  time-steps,  one  may  want  to  either  use  some  interpolation 
scheme,  or,  to  numerically  integrate  the  time-integrals. 

I'ne  author  'nas  developed  explicit  forms  ror  equation  (r.l)  by  assuming 
a  quadratic  interpolation  function  for  all  che  spatial  matrices.  However, 
ti:e  resulting  form  involves  a  large  number  of  matrices  which  m.ust  be 

numerically  evaluated  and  stored.  In  retrospect,  a  faster  method  appears 

to  be  using  program  and  input  logic  to  ring  certain  boundary  conditions 
■..‘hich  varv  ccnsiderablv  over  a  time-seen,  and  then  use  Gaussian  ?uadrature'' 


ccnsiceranl’. 


3. A.  Bole-)  and  J.F..  Weiner,  THEORY  GF  THERMAL  STRESSES,  Char.  L,  -John 
Wileii  ami  .Sans,  Ine.,  1960 

*  /I  -aussian  o-^der  W  will  exaozliJ  intearate  v  roliJKo^ia''  nf  mrat’v  I''-!. 
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to  numerically  integrate  the  corresponding  time-integral.  Hence,  this 
latter  approach  will  be  used  in  the  current  study. 


Because  of  the  desire  for  large  time-steps,  the  author  selected  an 
interpolation  function  for  a  quadratic  time-element,  i.e., 

][V(n)  ]  =  N^(2n)[V(2n)]  +  N^(2n+1)  [v(2n+l)  ]+  (2n+2)  [V(2n+2)  J  (F.5) 

where  N|.(2n)  =  -t(1-t)/2, 

Nj.(2n+1)  =  1-T 

N|.(2n+2)  =  t(1  +  t)/2, 

T  =  2(t-t  .  ,  )/At  , 
middle  n 


At  = 
n 


time-step  associated  with  the  nth  time-element. 


t 


middle 


t  +  At  /2 
n  n 


^V(2n)]  =  nodal  velocities  at  start  of  tiP.e-;-;teo , 
[v(2n-!-l)]  =  nodal  velocities  at  middle  of  time-step, 

and , 


[v(2n+2)]  =  nodal  velocities  at  the  end  of  the  time-step. 


If  we  now  assume  that  the  nth  weighting  function  is  non-zero  only  over 


the  nth  time-element,  then  t  =  At  and  the  first  term  in  equation  (F.4) 
can  be  rewritten  as. 


y"([KX]+[KYj+[KZj)[lV]WT^(t)dt  = 


n-1 


/wT  (t)  ([KX]+[KY]+[KZ] )I[U  ]+  L  At  J  (-  1^1-t) 
y  n  '  °  j=l  2  y  \  2 


d7[N(2j)] 


'n-l 


y"  (l-r“)dT[N(2j+l)]y  _r(l+T)dT[N(2j+2)]|  +y'  NT^(t)[v(n)]dt  Jdt. 


(F.6) 


-1 


■1  2 


n-1 


Substituting  for  the  quadratic  interpolation  functions  and  N  integrating 
equation  (F.6)  becomes  the  following  equation. 


V  "T 


[M(2n+1)]-  ^  /[M]WT  (T)fdT, 


At 


2  -1 


■M(2n+2)]=  2  /Vm]WT  (f)(l+2f)dT, 

,2-1  " 

At 

n 


[MF(2n)]=  -1  /[MF(f)]WT  (f)  {r ( W ) }dT  , 
. —  1  n 

At  -1 
n 

1  .  -7  _ 

[l^F(2n+l)]=  2  ;[MF(T  )]WT  (t.)  (1-T ')dT  , 

At  -1  ^ 

n 


CMF(2n+2)]=  1  /[t'lF(T)jWT  (f )  {fCl+T)  }d-r, 

—  1  n 

At  -1 
n 

1 

[kv1=;  WT  (i:)  ([KX]+[KY]+[KZ])df* 

-1  " 


2[U  ]  n-1  At. 

r  O  ,  ^  _ J_(iLV(2j)j 

_  3 

n  J*1  At 

n 


+4[V(2j+l)  ]+i[V(2j+2)  J)  } , 
3  3 


—  2  ^ 

[FT]=  t7  /[FT]WT  (T)dT, 

u  L  tl 

"-1 


and , 

[NL(k)]=  _J. 

At 

n 

In  the  above  equation,  rv(2n)]  corresponds  to  the  known  nodal  velocities 
at  the  start  of  the  time-step,  and  hence  there  are  twice  MT  number  of  un¬ 
knowns  and  equafions  represented  by  equation  (F,8).  This  is  the  disadvantage 
of  the  quadratic  time-element,  and  for  large  MT  considerably  more  computer 
time  is  required. 

2 

Equation  (F.8)  has  been  divided  by  the  ratio  (Atj^/2)2  to  eliminate  a 
similar  factor  which  appeared  as  the  coefficient  of  the  [K]  matrices.  .Also, 
the  spatial  matrices  have  been  retained  inside  the  time  integrals  so  that,  if 
necessary,  time  dependent  behavior  can  be  numericallv  integrated. 


[NL(k)]WT^(T)dT. 
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The  last  step  is  to  select  the  weighting  functions  WT  .  In  the 
Galerkin  approximation,  one  selects  the  weight  function  to%e  the  same 
function  as  used  for  the  element-interpolation.  This  implies  three 
functions  which  yield  three  matrix  equations  that  can  be  combined  into  a 
single  time-element.  The  single  time-element  is  then  combined  with 
successive  quadratic  time-elements  which  yields  the  recursion  relationships 
(see  Chap.  17  of  reference  8).  The  disadvantage  of  this  approach  is  Chat 
a  series  of  simultaneous  equations  must  be  solved  at  each  time-step.  To 
the  author's  knowledge  commercial  FE  programs  use  recursion  relationships 
which  assume  only  a  single  weighting  function  that  is  non-zero  over  a 
single-time  element.  Furthermore,  most  programs  use  linear  time-elements 
which  yield  the  familiar  forward  difference,  central  difference  and  back¬ 
ward  difference  type  of  equations. 

If  the  weighting  function  weights  over  more  than  one  time-step,  and 
a  linear  time-element  is  assumed,  then  higher  level  schemes,  such  as 
HoufaoLt's  Method,*  are  obtained.  The  advantage  of  Che  higher  level  scheme 
over  the  quadratic  time-element  is  that  there  are  only  MT  unknowns  per  time- 
step.  However,  it  appears  chat  the  number  of  stored  matrices  for  the  three- 
level  scheme  is  about  Che  same  as  that  for  the  quadratic  time-element. 
Therefore,  computational  speed  involves  a  trade-off  between  more  storage  and 
CPU  time  per  time-step  for  the  quadratic  time-element  versus  fewer  time- 
steps  required  to  solve  a  problem. 

One  should  note  that  the  usual  displacement  formulation  of  the  linear 
momentum  equation  involves  a  second-order  equation  in  time,  and  a  quadratic 
time-element  is  used  in  many  current  FE  program®.  The  corresponding 
recursion  relationship  is  a  three-level  scheme^  with  a  variety  of  forms 
similar  to  Newmark's  general  algorithm.^^  Some  FE  programs  use  a  four- 
level  scheme,  where  there  is  still  only  MT  unknowns,  and  special  "starting" 
algorithms  must  be  used  to  determine  the  unknowns  at  the  first  three  time- 
nodes. 

In  order  to  perform  numerical  experiments,  the  author  assumed  that  the 
spatial  matrices  of  Che  left-hand-side  of  equation(F.8)  are  independent 
of  time  within  a  time-step,  and  equation  (F.8)  was  rewritten  in  the 
following  form, 

CAj[V(2n+2)]+[B][V(2n+l)]+[Cl[V(2n)]=[CF]  ,  n*0,l,...NT  (F.9) 


* 

See  Swanson  Analysis  Systems,  Ina,  Theoretical  Manual  for  application  of 
this  method  to  a  vanity  of  FE  programs,  2nd  Ed,  1982 

JO 

'ILM.  ^lewmark,  "A  Method  for  Computation  of  Structural  Dynamics" ,  Free. 
Am.  Soc.  Civil  Eng.,  35,  EMS,  pp.  67-94,  2959. 
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where  the  [a],[b],[c]  and  [CF]  matrices  take  different  forms  depending 
upon  the  weighting  function  assumption. 


When  a  single  weighting  function  is  assumed  over  the  time-step, 
equation  (F.9)  becomes  the  usual  three-level  scheme  with  the  following 
change  of  the  indices;  2h+'2-^n+l ,  2n+l->-n,  and  2n->-n-l ,  with  n=l,  2,..., NT 
Table  F-1  lists  six  weighting  functions  selected *by  the  author  and  the 
corresponding  matrix  elements.  Note  that  the  [MF]  terms  have  been 
neglected , Ck]=[KX]+[KYJ+[ KZ] ,  and  the  [f]  matrix  reoresents  the  right- 
hand-side  of  equation  (F.3).  All  matrices  are  evaluated  at  the  middle 
of  the  time-step  At. 


When  three  weighting  functions  are  assumed  over  the  time-step 
the  following  form  of  the  weighted  residuals  is  obtained  for  a  single 
quadratic  time-element, 


^11  ^12  ^13 

'  V(2n) 

•KinF(2n) 

^21  ^22  ^23 

V(2n+1) 

= 

K,^F(2n+l) 

1/  If  V 

31  32  33 

V(2n+2) 

K2^F(2n+2) 

- 

(F.IO 


where  the  K..  values  are  dependent  upon  the  assumed  weighting  functions. 
In  this  set^if  equations,  the  V(2n)  are  known  (either  the  initial  values, 
or,  from  the  previous  time-step).  Therefore,  equation  (F.IO)  can  be  re¬ 
written  as 


^22  ^23 

’V(2n+1)" 

F(2n+l)-K2^V(2n)' 

.^32  S3_ 

V(2n+2) 

K.  F(2n+2)-K.,V(2n) 
u  3n  31 

2MTxl 


As  discussed  earlier,  this  larger  system  of  equations  is  justified 
only  if  either  a  better  approximation  is  obtained  for  the  same  CPU  time, 
or  the  same  approximation  is  obtained  with  less  CPU  time  than  the  single 
weighting  function  approach. 

In  order  to  evaluate  the  above  K. .  coefficients,  the  author  assumed 
the  following  two  setsrA.  A  true  Gale^kin  approximation  using  the  second 
through  the  fourth  weighting  functions  given  in  Table  F-1,  and  B.  three 
linear  functions  associated  with  the  last  two  weighting  functions  given 
in  Table  F-l_^  The  resulting  K..  coefficients  are  given  in  Table  F-2 . 
Again,  the  [MF]  terms  have  been"^ neglected . 


* 

These 


are  the  sane  as  those  used  in  reference  ~5  sc 
Trade  itith  the  standard  disrlacer.ent  approach. 


Table  F-1  Possible  Weighting  Functions  and  the 
Corresponding  MTxMT  Matrices 


Table  F-2  Coefficients  for  Equation  (F.ll) 
Associated  with  Two  Sets  of 
Weighting  Functions 


Coefficient 


Set  A 


Set  B 


ll[K(2n)]-  2[M(^)] 


3[K(2n)]+[M(2n)] 


2[K(2n+l)] 

9 


[K(2n+1)] 

6 


[K(2n+2)]  +  2fM(2n+2)1 


-rK(2n+2)]+[M(2n+2)] 


K  ll[K(2n)MM(2n)] 

3,^2 


43[K(2n)]+[M(2n)] 


Koo  ^^[K(2n+l)]-4[M(2n+l)] 

45  2 

3At 


3[K(2n-fl)]-4[M(2n+I)] 


CK(2n+2)]+  [M(2n+2)] 


13[K(2n+2)]+7[M(2n+2)] 


2.  Fluid-Type  Behavior 

As  discussed  in  Appendix  E,  Che  most  complicated  matrix  involves  the 
strain-rate  matrix  of  equation  (E.13).  Inspection  of  this  equation  shows 
that  each  array-element  involves  the  velocity  matrix  [v].  Furthermore, 
since  each  matrix  product  within  the  (UVW...)  type  of  array-element 
corresponds  to  a  scalar,  the  (UXV)  (UZIV)  type  of  items  can  be  rewritten 
as  (UZIV)  (UXV),  etc.  Therefore,  each  array-element  can  be  considered 
as  post  multiplied  by  [v]  and  equation  (E.13)  can  be  rewritten  as 

[NLU]  =  [UXYZ]  [V].  (F.12) 

6xMT 

where  [UXYZ(l)]  =  2  { (UXIV)  [UX]  +  (VXIV)  [VX]  +  (IJXIV)  [WX]  }, 

IxMT 

[UXYZ(2)]  =  2  {(U’YIV)  [UY]  +  (VYIV)  [VY]  +  (WYIV)  [WY]  }, 

[UXYZ(3)]  =  2  {(UZIV)  [UZ]  +  (VZIV)  [VZ]  +  (WZIV)  [WZ]  }, 

[UXYZ(4)]  =  (UXIV)[UY]+(LTIV)[U'X] 

+(VXIV) [VY]+(VYIV)[VX] 

-(-(raiv)  ]+(WYIV)  [ra] , 

LUXYZ(5) ]-(UXIV)[UZ]+(UZIV)[UX] 

+(VXIV)[VZ]+(VZIV)rvX] 

4-(WXIV)  [WZ  ]+(WZIV)  [wx  J  , 

EUXYZ(6)]=(U^IV)[UZ]+(UZIV)[UY] 

+(VYIV)[V2]+(VZIV)[VY] 

+(WYIV) [WZ]+(WZIV) [WY] . 

Multiplying  equation  (E.16)  by  the  nth  weighting  function,  substituting 
equations  (F.3),  (F.5)  and  (F.12),  integrating  over  time  and  rearranging 
terms  yields  Che  following  equation, 

([i^(2n+2)]+[G(2n+2)]+[M(2n+2)]-[MF(2n+2)]+CK^(2n+2)]  ) 

[V(2n+2)  ]+([^(2n+L)  ]+CG(2n+l)  ]+[M(2n+l)  ]-[MF(2n+l)  ] 

+[KEN(2n+I)  [)[V(2n+l)  ]+(L^(2n)  ]+[G(2n)  ]+[M(2n)  ]  (F.  13) 

-[MF(2n)-]+[i^(2n)  J)  [V(2n)  ]=-[GV]-[FB]+[FT] 

_  _  _ 20  _ 

+[ PE  ]+[ GXYZ  >r  ET AXYZ  ] - ![ NLF  (  k)  ] 

k=L 

k  4=  4,5,6 
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where  [KE(2n)]  =  1  /  ([KEX]  +  [KEY]  +  [KEZ])  WT  (t){-t(1-t) }dT, 

aT  -1 

n 


[KE(2n+l)]  =  2  /  ([KEX]  +  [KEY]  +  [KEZ])  WT  (T)(l-T^)dT. 

I -  1  ri 

At  -1 
n 


_  1 

[KE(2n+2)]  =  /([KEX]  +  [KEY]  +  [KEZ])  WT  (T)T(l+T)dT, 

At  -1 

n 

_  1 

[KEN(2n)]  =  1  /([KENX]  +  [KENY]  +  [KENZ])  WT  (t ) { -T ( 1-t) Idi , 

At  -1 

n 


[KEN(2n+l)]  =  2  /  ([KENX]  +  [KENY]  +  [KENZ])  WT  (T)(l-T^)dT, 

At  -1  " 

n 

_  1 

[KEN(2n+2)]  =  /  ([KENX]  +  [KENY]  +  [KENZ])  WT  (T)T(l+T)dT, 


[G(2n)  J=/[GXYZ](_5.-il+ii)WT  (-c'/dT, 
-1  12  4  6  " 


—  ^  «3  -  - 

[G(2n+l)j=/[GXYZ](2+T-T  )WT  (T)dT, 

-1  3  3” 

[G(2n+2)]=;rGXYZ](2i+ii  +t2)WT  (T)d7, 

-1  12  4  6  ” 


_  1 

[GVl=.rWT  (T)[GXYZ]dT*{2[U  ] 
,  n  o 


Z  j5^(l[V(2j)]+4[V(2j  +  l)] 

i=l  At  3  3 

.  n 

+i[V(2j+2) ]) } 

3 


_  1 

[PE]«  2  ;{[PEXMPEY]+[PEZ]}WT^(T)dT. 

At  -1 
n 

[cm  J=_2  ;  {[CX]''^[GX]  +  [CY]^[GY]  +  [CZ]^[GZ] }  [e^]dVol  WT^(T)d7, 

At  -1 
n 

_  1  _ 

[ETAXYZ]=_2  /[ETAXYZ]WT^(T)dT, 

At  -1 
n 

_  1 

[NLF(k)]=^  /  [NLF(k)]WT  (T)dT, 

At  -1 
n 

[  KEX  ] =/[  CX  ]"“■  [  ETAX  ]  [  B  ]  d  Vo  1 , 

V(j) 

[KEY]=/[CY]'^[ETAYl[B]dVol, 

V(j) 


[KEZ]=/[CZ]'^[ETAZl[B]dVol, 

V(j) 


[KENX]=/[CX]^([I]+Ca])[ETAX][UXYZ]  +  [A][ETAX][B]}dVol, 
V(j) 

[KENY]»;[CY]^{([I]+[a])[ETAY][UXYZ]  +  [A][ETAY][B]}dVol, 
V(j) 


[KENZ]=/CCZ]'^{([i]+[a])[ETAZ][UXYZ]  +  [A]  [ETAZ]  [B]  }dVol, 
*  V(j) 

[PEX’=/(p+p  )CUX]^  dVol, 

V(j)  ® 

[PEY]=/(p+p  )[UY]'^  dVol, 

V(j)  ® 

[PEZ"=/(p+p  )[UZ]^  dVol, 

V(j)  ® 


[Gm>/((CX]^[GX]  +  [CY]'^[GY]  +  tCZ]’^[GZ])[B]dVol, 
V(j) 


m 


>.y. 
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[  ETAXYZ  ]=/ (  [  CX]'^[  ETAX  ]+[  CY  ]'^[  ETAY  ] 

V(j)  .  5. 

+[CZ]  [ETA2])[e  ]dVol. 

Note  chat  the  nonlinear  terms  NLF(k),  k  =4 . 5 . 6 . 2 1 .  22  and  23, 
in  equation  (E.16)  have  been  rearranged  into  the  [KEN]  terms  of  the  above 
equation.  Also,  equation  (F.13)  has  exactly  the  same  form  as  equation 
(F.8),  except  that  the  velocity  coef f icient-matrices[KE]  and  [KEN]  depend 
upon  the  unknovm  velocity.  Hence,  even  if  all  Che  [NLF]matrices  are 
neglected,  equation  (F.13)  can  only  be  solved  via  an  iterative  process. 

_  For  the  special  case  of  no  strain-hardening  i.e.,  [GXYZ]eO, 
the  [g]  and  [GV]  matrices  are  identically  zero  and  equation  (F.13)  contains 
a  common  At  term  which  can  be  canceled  out.  The  resulting  equation 
corresponds  to  the  familar  first-order  type  of  recursive  equation. 

Since  the  form  of  equation  (F.13)  is  similar  to  equation  (F.8), 
the  explicit  equations  for  the  same  assumptions  and  specific  weighting 
function,  or  set  of  weighting  functions,  will  be  the  same.  Therefore,  the 
coefficient  matrices  given  in  Tables  F-1  and_F-2  also  apply  to  equation 
(F. 13)  if  the  [k]  matrix  is  replaced  by  the[GXYZ  ]  matrix,  and  the  [  KE j, 

[ken  ],  and [  MF  ]  matrices  are  neglected. 

For  the  more  general  problems,  the  [K]  matrix  of  Tables  F-1  and 
F-2  must  be  replaced  by  three  matrices  corresponding  to  [Gl,[KE]  and  [KEN]. 
Using  the  same  functions  as  previously  for  the  elastic  type  of  behavior, 
Table  F-3  and  F-4  list  the  corresponding  [Aj,[B]  and  [C]  matrices  for  usage 
in  equation  (F.9).  Similarly,  Tables  F-5  and  F-6  list  the  K. .coefficients 
for  usage  in  equation  (F.ll)  when  applied  to  fluid-type  behavior. 


Table  F-3  Possible  Weighting  Functions  and  the 
Corresponding  MTxMT  Matrices  for 
Fluid-Type  Behavior  and  Equation  (F.9) 

Function:  Constant  =  1 

[a]=  1  ([ke]+[ken])+[m] 

At^ 

lB]=  2  ([KE]+[KEN])4-[GXY2] 

3At  3 

[c]=  1  ([ke]+[ken])+[gxyz]-[m] 

6At  6  .J2. 

At 


Function:  -t(1-t)/2 

Ca]=  -1  ([ke]+[ken])-[gxyz]-[m] 

5At  60  .2 

At 


[b]=  2  ([KE]+[KENj)+13[GXYZ]+  4  [m] 

5At  30  ^^2 


[C]*  4  ([KE]+[KEN])+[GXY2]-  3  [M] 

A  r  19  2 


Function:  (1-t  ) 

[a]"=  _ ([I^1+CKEN])-[GXYZ]+[Mj 

lOAt  60  ^^.2 

[B]=  4  ([KE]+[KEN])+[GXYZ] 

5At  3 

[c]=_l _ ([KE]+lKEN])+JJ[GXYZ]-l^ 

lOAt  60  ^2 


Table  F-4  Possible  Weighting  Functions  and  the 

corresponding  MTxMT  Matrices  for  Fluid 
Type  Behavior  and  Equation  (F.9) 


Function:!  (1+t)/2 

[a]=  4  ([KE]-f[KEN])+[GXYZ]+ 3[M] 

5At  12  ^^2 


Cb]=  2  ([KEMKEN])+n[GXYZ]-  4  [m] 
5At  15 

At 


[c]=  -1  ([KE]+[KEN])+11[GXYZ]+.IM]. 

5At  60  .,2 

At 


Function:  -t,  t  ,<t<t  ;t,  t  <t<t 


[A]=  1  (rKEl4fKENl)4f  GXYZI-H'MI 

4At  48  2 


[B]=  1  ([KE]+[KEN])+[GXYZ] 

2At  3 


[C]=  J[_([KE]+[KEN])+_2[GXYZ]-JM]_ 
4At  48  ^j.2 


Function:  l+x,  t  ,<t<t  ;  1-t,  t  <t<t 

n-1-  -  n  n-  - 

[A] =  _1 _ ([ke]+[ken])-[gxyz]+[m] 

12At  96  ^,2 

At 

[B] =  5  ([ke]+[ken])+[gxyz] 

6At  3 


[ C ] -  1  ( [ KE ]+[KEN ] ) +_3l GXYZ ] - [  MJ 


12At 


16 


At* 
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Table  F-5  Coefficients  for  Equation  (F.ll) 
Associated  with  Set  A  of  the 
Weighting  Functions 


Coefficient 


II  [G(2n)  ]+  J^([KE(2n)  ]+[KEN(2n)  ])-  2[M(2n) 
90  l5At  2 


i[G(2n+l)]+  _3_  ([KE(2n+l)]+[KEN(2n+l)  ]) 
9  15At 

n 


'23  -[G(2nj^]+  l_([KE(2n+2)]+rKEN(2n+2)])  +  2[M(2n+2)] 


90  15Ac 


K  n_[G(2n)]  -_1 _ (rKE(2n)  MKEN(2n)  ]-f[M(2n)  ] 

180  L5At^  3^^  2 

n 

ll[G(2n+l)]+  2  ([KE(2n4-l)]+[KEN(2n+l)])-4[M(2n+l)] 

45  15At  2 

n  3At 

n 


K,,  CG(2n+2)]+_4___([KE(2n+2)]+CKEN(2n+2)])+  [M(2n+2)] 

-16-  15Ar„  2 


K,  2/(36t  ) 

2n  n 


K-  l/(3At  ) 

in  n 


G(2n) ,  G(2n+1),  G(2n+2)  =  [GXYZ]  evaluated  at  times 

•^Zn,  hn+l,  "2n-2  respectively. 


Table  F-6  Coefficients  for  Equation  (F.ll) 
Associated  with  Set  B  of  the 
Weighting  Functions 


Coefficient 


21 


3[G(2n)]+  1 
32  24At 


(rKE(2n)  MKEN(2n)  ])-f-[M(2n)  ] 

2At  ^ 


K,,  [G(2n+1)  ]+_5_([KE(2n+.)  ]+[KEN(2n+l)  ]) 

6  12  t 

n 


K  -[G(2n+2)]+  1  ([KE(2n+2)  MKEN(2n+2)  l)+[M(2n+2)  ] 

192  24At^  2,^  2 


K  «  G(2n)  -J. _ (rKE(2n)  MKEN(2n)  ])+[M(2n)  ] 

480  24At  .,^2 

n  6  At 

n 

K  3 [ G ( 2n+l ) ]+  1  ([KE ( 2n+l ) ]+[ KEN ( 2n+l ) ] ) -  4  [M ( 2n+l ) ] 

in  AAt-  -  2 


K 


33 


13  [G(2n+2)]+  7 
480  24At 

n 


( [KE ( 2n+2 )  ]+[ KEN ( 2n+2 )  ] ) -f  7  [M(2n+2) 

6At  ^ 
n 


K_  l/(2At  ) 

zn  n 

K_  l/(2At  ) 

3n  n 
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The  discretized  linear  momentum  equation  contains  coefficients  and 
a  pressure  term  which  depend  upon  the  temperature.  The  temperature 
explicitly  appears  in  the  energy  equation  (D.4),  and  therefore,  the 
momentum  and  energy  equations  are  coupled  together.  The  fully  discretized 
velocity  can  be  utilized  directly  in  equation  (D.4),  and  the  temperature 
is  spatially  discretized  using  the  following  approximation, 

NR  NPT(j) 

T(x,y,z,t)  s  E  z  NTE  (x,y,z)  T.(t)  (G.l) 

j=l  i»l 

NR 

3  Z  [NTE]  [T(t)] 

j  =  l 

where  NTE  =  a  trial  function  in  terms  of  the  spatial  variables, 

T^(t)=  the  unknown  value  of  the  temperature  at  the  ith  node, 

NPT(j)=  number  of  temperature-nodal  points  in  the  jth  region. 

Equation  (23)  is  spatially  discretized  by  substituting  equation  (G.l) 
to  yield 

5  (thermal  part  of  EOS)  =  +11^  [NTE][T] 

3  t  dT  3T 

+  B.  +  ^  [T]’^  [NTE]”^  [NTE]  [T]  (G.2) 

3T 

+  2B^  [NTE][T]}[NTE][f]  +  (_3_S^  [NTE]'^)  [NTE]  [T] . 

3 1  syfj  3  t 


The  spatial  gradients  of  the  temperature  are  discretized  into  the 
following  forms 

K  3T  n  +  K  3T  n  +  K  3T  n  =  (K  n  [TX] 


+  K  n  [TY]  +  K  n  [TZ])  [T] , 
yy  y  zz  z 

VW'(k  3T  +  k  3T  +  K  3T)  =  [BT]^  [KT]  [BT]  [T] , 
XX—  yyT~  zz— 

3x  3y  3z 

where  [TX]  =  [^  NTE^  ...  _3_  NTE^j^p]  , 

3x  3x 


(G.3) 


(G.4) 


[TY]  =  [_d_  NTE.  .  .  .-  _3_  NTE^j^p] , 
3y  3y 


[TZ]  =  [  9  NTE,  ...  9  NTE_], 


[BT]  = 


[TX]  [TY]  [TZ]]^j^p^3 


The  last  term  that  must  be  discretized  is  the  product  involving  the 
strain-rate  and  stress.  Using  previously  defined  matrices,  the  strain-rate 
matrix  becomes, 

=  [V]^  [B]^  +  [V]”^  [UXYZ]'^.  (G.5) 

The  stress  matrix  depends  upon  the  condition  of  the  material.  For 
elastic  behavior,  the  stress  matrix  can  be  written  as, 

L^iJ  =  [e]  -  [p].  (G.6) 

■‘NFGxI  -'nfgxnfg 


Therefore,  the  matrix  product  for  elastic  behavior  is 

=  [V]'^  ([B]^  +  [yXYZ]”^  (G.7) 

{[C..]  ([B]  [IVl  +  (NLUl)  -  [p]}. 


For  fluid-type  of  behavior,  equations  (15)  through  (17)  are  used  to 
yield  the  following  stress  matrix  for  loading  behavior, 

[a.^]  =  [G]  ([c]  -  [e^])  +  [ETA]  ([e]  -  [t^])  -  [p^l ,  (G.8) 

where  [G]  =  G 

[0]  [I] 


*The  previously  defined  matrices  [D^],  i=x,y,z,  are  for  the  special 
case  of  an  isotropic  material. 
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[ETA]  = 


2[I]  [0] 
[0]  [I] 


[p^]  *  [p^p^p^OOO]  , 


Px  “  P  +  Pe> 


when  NFG  equals  six.  The  corresponding  matrix  product  becomes, 

[e]’^  [a..]  =  [V]'^([B]'^  +  [UXYZ]'^)|[G]([B]  [IV]  +  [NLU]  -  [e^]) 

+  [ETA]{([B]  +  [UXYZ])  [V]  -  [c^]}  -  [p^]}.  (G.9) 

If  we  now  assume  that  the  weighted  momentum  equations  are  zero  for  the 
discretized  values  of  temperature  and  velocity,  and  the  weighting  function 
equals  the  trial  function,  energy  equation  (D.4)  can  be  written  in  the 
following  form  for  each  region, 

{  fo  (a,(T)  +  a,(T)  [NTE]  [T]  +  [T]^  [NTE]"^  [NTE]  [T]) 


[NTE]^  [NTE]dVol}[T]  +  +  a^(T)  [T]^ 

V(j) 

[NTE])  [NTE]^  [NTEjdVol  +  j['S,T]  [KT]  [BT]dVol  +  jh  [NTE]^  [NTE]dS 


S,(j) 


+  /a(T,T  )  [NTE]^  [NTE]dS}  [T]  =  /[NTE]^  qdS 

(u)  "  s^U) 

+  y’p[NTE]^  QdVol  +  yE  T^  [NTE]^dS 
V(j)  S^(j) 

-  j  [NTE]^  (Internal  energy)  (p[N]^  [V]’n)dS 

Sp(t) 

f  T  5 

+/  [NTE]  a(T,T  )  T  dS  +  E  [FV(k)], 

I  00  00 

•c  l■\ 

S^(j) 


(G.IO) 
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where  a(T,T  )  = 
a^(T)  = 

32(1)  = 

33(1) 

3^(T)  = 

a  = 

e 

°SB 

T 

00 

h 

S 

q 

s^Cj)  = 

SfaCj)  = 

s,a)  = 

[FV(1)]  = 
[FV(2)]  = 
[FV(3)]  = 
[FV(4)]  = 
[FV(5)]  » 


a^FfOgg  ([T]^  [NTE]”^  [NTE]  [T]  +  )*  ([NTE]  [T]  +  T J , 

^  +  B  , 
dT  ^ 

2B  +  ^  , 

3T 

Ml  1^* 

3vT^  3t 

M2 

3t 

effective  emissivity, 
form  factor, 

Stephen-Boltzmann  constant, 
far-field  environmental  temperature, 
effective  film-coefficient, 

t 

s  +  s,  +  s  , 

a  b  c 

boundary  of  jth  region  which  has  a  conductive  flux  q, 

boundary  of  jth  region  which  has  a  convective  flux, 

boundary  of  jth  region  which  has  a  radiation-flux, 


-/pInte]”^ 

v(j) 

-;p[NTE]'^ 

V(j) 

/  [NTE]^ 
S(j) 

/  [NTE]^ 
S(  j ) 

/p[NTE]^ 

V(j) 


[V]"'^  [N]”^  [FgldVol, 

[V]^  [N]^  [VjdVol, 

[V]^  [N]^  [n]  [ST]dS, 

[V]''^  [N]^  [A]  [n]  [ST]dS, 
(fp-f^)  [e]^  [o^^JdVol. 


The  flux  q  is  either  specified  on  an  exterior  boundary  of  the  body,  or 
evaluated  from  an  adjacent  region  using  equation  (G.3).  The  product 

[e]  [a..]  is  given  by  either  equation  (G.7)  or  (G.9)  depending  upon  the 

materiai^ condition .  The  velocity  variable  is  implicitly  contained  within 
the  coefficients  of  the  left-hand-side  of  equation  (G.IO),  and  explicitly 
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contained  within  the  right-hand-side  matrices  [FV(k)],  k=l,...4.  The 
temperature  variable  is  implicitly  contained  within  the  [FV(5)]  matrix. 


Note  that  the  usual  quasi-linearization*  of  the  radiation  term  has  been 
used  in  equation  (G.IO).  However,  even  without  the  radiation  term,  equation 
(G.IO)  is  an  extremely  non-linear  ordinary  differential  equation  and  must  be 
solved  using  some  sort  of  iterative  procedure.  When  the  and  B2 

contributions  are  ignored,  then  aj^(T)  =  c^,  a^CT)  =  a^CT)  =  a^(T)  E  0  and, 

if  there  is  no  radiation,  the  left-hand-side  of  equation  (G.IG)  takes  the 
more  familiar  form  of  pc^[T]  +  [KT][T]. 

Excluding  the  implicit  time-dependency  of  the  a^(T)  coefficients,  the 
left-hand-side  of  equation  (G.IO)  can  be  rewritten  in  the  following  form  of 
temperature-dependent  spatial  integrals  times  the  time-dependent  temperature 

L.H.S.  -1-  ([CO]  +  [Cl]  +  [C2])  [T]  +  ([KO]  +  [Kl]  +  [K2]  +  [K3] 

+  [K4])  [T],  (G.li) 

where 

[CO]  =  /  pa  (T)  [CT]dVol, 

V(j) 

[Cl]  *  f  pa,(T)  b,(T)  [CT]dVol, 

V(j) 

[C2]  *  /  p  3B,  bf(T)  [CT]dVol, 

V(j)  oT^ 

[KO]  =  /  [BT]  [KT]  [BT]dVol, 

V(j) 

[Kl]  =  /  h  [CT]dS, 

S,(j) 

[K2]  =  /  a(T,T  )  [CT]dS, 

00 

S^(j) 

[K3]  =  ;  pa^d)  [CT]dVol, 

V(j) 

[K4]  =  /  pa  (T)  b,(T)  [CT]dVol, 

V(j) 

[CT]  =  [NTE]^  [NTE], 

^  NPT(j) 

b,(T)  =  [NTE]  [T]  =  E  NTE  T.(t). 

^  i=l 

*(t'^-T  ‘*)  =  (T^  +  T  (T  +  T  )  (T  -  T  )  -<■  a(T,  T  )  (T  -  T  ). 
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Note  that  a  mild  nonlinearity  is  introduced  by  the  [CO]  and  [K3]  terms 
whereas,  stronger  nonlinearities  are  introduced  respectively  by  the  [Cl], 
[K4],  [C2]  and  [K2]  terms.  The  volume  and  surface  integrals  cannot  be 
numerically  evaluated  until  the  temperature  is  known,  and  hence  these 
temperature-dependent  spatial  integrals  can  only  be  determined  in  some 
iterative  manner. 

In  order  to  discretize  equation  (G.lO)with  respect  to  time,  equation 
(G.lO)is  multiplied  by  a  piecewise,  non-zero,  time-dependent  weighting 
function  WTE(t)  and  a  trial  function  is  introduced  to  approximate  the 
temporal  temperature  behavior,  i.e., 

NT 

WTE(t)  =  Z  WTE  (t),  (G.12) 

m=0  “ 


[T(t)]  3  Z  N  [T(m)], 

m=0 


(G.13) 


where  [T(m)]  represents  the  TNP  nodal  temperatures  at  time  t^,  and  N^^  is 

the  time-dependent  interpolation  function  for  the  temperature.  Integrating 
over  the  entire  time  interval  of  interest,  and  using  equation  (G.ll),  the 
energy  equation  becomes  the  following  set  of  nonlinear  algebraic  equations 
for  each  time-step, 

{  f  J_  ([CO]  +  [Cl]  +  [C2])  WTE^  (t)  [N^gld? 

J  At 
c  m 
s 


+  /  ([KO]  +  [Kl]  +  [K2]  +  [K3]  +  [K4])  WTE^(t)  [N^.^]  dx} 


(G.14) 


[T(m) 


]  =  /  {[q]  +  [Q]  + 


[h]  -  [IE]  +  [RD]  +  Z  [FV(k)]}WTE  (T)dT, 

k=l 


where  [q] 


I  [NTE]  q  dS, 
S^(j) 


[Q]  =  /p[NTE]  Q  dVol, 

V(j) 

[h]  =  /h  T  [NTE]^  dS, 

[IE]  =  /  [NTE]^  (Internal  energy)  (p[N]^  [V]-n)  dS, 

Sp(t) 

■RD]  -  /  [NTE]"^  a(T,TJ  T^dS. 

S^(j) 


_  Note  that_the  [FV]  terms  include  the  time-dependent  function 
[V]  s  [V(n)],  which  is  explicitly  given  by  equation  (F.5).  Also, 

the  time-step  At^  for  the  thermal  approximation  could  be  larger  than  the  At^ 

used  in  the  velocity  approximation,  in  which  case  the  [V]  term  would  involve 
a  sum  of  interpolation  terms. 

If  a  quadratic  time-element  is  assumed  for  the  temperature  behavior, 
then  the  algebraic  equations  take  the  form  of  equations  (F.9)  or  (F.ll). 
Similarly,  the  coefficients  are  given  by  Tables  F-3  through  F-6,  with  the 
[GXYZ]  and  [G]  matrices  omitted,  all  coefficients  multiplied  by  At  ,  and  the 
following  substitutions  made:  ([KE]  +  [KEN] ) ->  [KT] ,  and  [M]  -*■  [CT]^  where 

m  =  /([KO]  +  [Kl]  +  [K2]  +  [K3]  +  [K4])  WTE^  (7)  [N^^]  d7, 

^s 

[CT]  E  /([CO]  +  [Cl]  +  [C2])  WTEjjj  (7)  [N^^]  d7. 

^■s 

For  ballistic  penetration  problems,  the  temperature  change  between 
time-steps  might  be  realistically_modeled  using  a  linear  variation.  If 
a  single  weighting  function  WTE^(t)  is  assumed  over  the  time  interval,  then 

equation  (G.14)  becomes  the  classical  "implicitness  0"  equation  (0<9£l)  for 
a  first  order  equation,  i.e.. 


(G.15) 


+  (1-9)  t>lT(T„^)l  >T„  ■ 


where 

and  which  can  be  solved  for  using  a  linear  iteration. 

4 

Rather  than  the  above  linear  iteration.  Baker  uses  a  Newton-Raphson 
iteration  scheme  of  the  following  form. 


where 


t-KVl”'’ 

-  iVi 

tVii' 

tVi>’' 

(G.16) 


tT„l)  +e  (1-9) 


[KT  ]  [T  ]  +  9  [F„_^,]  +  (1-0)  [F„], 


3[H]/3[T]. 

[^]  +G(^]  -k-f  1  S[^]  +  +  9[F] . 

\Atn  9[T]  3[T]y  "  ^  3[T] 

During  the  second  year  of  this  research,  several  first  order, 
non-linear,  equations  will  be  evaluated  using  a  variety  of  the  above 
techniques,  and  several  of  the  most  promising  techniques  will  be  selected. 
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